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SUMMARY 


This  is  the  Final  Report  of  the  research  project  ''Advanced  Quantitative  Robust  Control 
Engineering:  New  Solutions  for  Automatic  Loop-Shaping  for  SISO  and  MIMO  Systems.  Part  I: 
SISO  Systems.^'  (FOARD  Award  #  FA  8655-08-1-3027;  Grant  #  083027). 

The  outcome  is  a  new  automatic  loop-shaping  methodology  to  design  QFT  controllers  for  Single- 
Input-Single-Output  (SISO)  systems,  which  uses  two  techniques:  Fvolutionary  Algorithms  and 
Genetic  Algorithms.  The  methodology  is  based  on  the  classical  theorems  that  define  the  control 
limitations  on  performance  in  SISO  systems,  on  advanced  techniques  to  search  for  optimum 
solutions  subject  to  restrictions  defined  as  quadratic  inequalities,  and  on  previous  experience 
designing  QFT  compensators  on  the  Nichols  Chart.  Using  this  procedure,  it  is  possible  to  avoid 
classical  restrictive  assumptions,  such  as  rectangular  boundary  for  the  templates,  pre-specified 
types  of  bounds,  some  controller  parameter  fixed,  available  initial  controller  to  start  from,  etc.,  as 
well  as  restrictive  requirements  of  other  well-known  optimisation  techniques,  such  as  continuity, 
convexity,  or  differentiability.  The  new  approach  does  not  need  such  assumptions,  overcoming 
the  disadvantages  of  the  previous  automatic  procedures,  which  results  in  very  conservative 
designs,  and  is  easily  implemented  on  existing  CAD  tools.  In  conclusion,  it  contributes  to 
improve  QFT  controller  design,  alleviating  much  of  the  manual  work  required  at  present. 
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1.  INTRODUCTION 

The  Quantitative  Feedback  Theory  (QFT)  method  ([27-30]  and  [40])  offers  a  direct,  frequency- 
domain  based  design  approach  for  tackling  feedback  control  problems  with  robust  performance 
objectives.  In  this  approach,  the  plant  dynamics  may  be  described  by  frequency  response  data,  or 
by  a  transfer  function  with  mixed  (parametric  and  non-parametric)  uncertainty  models.  One 
feature  that  distinguishes  QFT  from  other  frequency-domain  methods  is  its  ability  to  deal  directly 
with  uncertainty  models  and  robust  performance  criteria.  This  is  achieved  by  translating  robust 
performance  specifications  and  uncertainty  models  into  so-called  QFT  bounds.  These  bounds, 
typically  displayed  on  a  Nichols  chart-like  plot,  then  serve  as  a  guide  for  shaping  the  nominal 
loop  transfer  function  which  involves  the  manipulation  of  gain,  poles  and  zeros.  This  design 
process  is  executed  efficiently  using  computer-aided  design  software,  such  as  the  QFT  Control 
Design  MATLAB  Toolbox  ([3]),  the  AFIT  CAD  tool  ([31],  [35]  and  [36]),  Qsyn  ([21])  and  [30], 
and  is  effective  for  simple  problems.  Nevertheless,  QFT  designers  are  often  challenged  by  such 
control  problems  due  to  a  lack  of  loop-shaping  experience,  and  could  benefit  from  an  algorithm 
that  automatically  provides  a  first  cut  solution  to  the  loop-shaping  problem.  In  addition,  an 
automatic  loop-shaping  facility  would  enhance  the  capabilities  of  the  expert  QFT  designer. 
Automatic  loop-shaping  algorithms  have  been  proposed  over  the  past  twenty  years  and  this  work 
reports  on  a  new  version. 


2.  AUTOMATIC  LOOP-SHAPING  METHODOLOGY 


2.1.  OVERVIEW 

Quantitative  Feedback  Theory  ([27-30]  and  [40])  is  a  robust  frequency  domain  control  design 
methodology  which  has  been  successfully  applied  in  practical  problems  from  different  domains 
[12-21].  One  of  the  key  design  steps  is  loop  shaping  of  the  open  loop  gain  function  to  a  set  of 
restrictions  given  by  the  design  specifications  and  uncertain  model  of  the  plant.  Although  this 
step  has  been  traditionally  performed  by  hand,  the  use  of  CAD  tools  ([3],  [30],  [31],  [35],  [36]) 
has  made  the  manual  loop  shaping  much  simpler.  However,  the  problem  of  automatic  loop 
shaping  is  of  enormous  interest  in  practice,  since  the  manual  loop  shaping  can  be  hard  for  the  non 
experienced  engineer,  and  thus  it  has  received  a  considerable  attention,  especially  in  the  last  two 
decades. 

The  concept  of  controller  design  automation  in  QFT  was  introduced  by  Gera  and  Horowitz  [22] 
who  proposed  an  iterative  procedure  based  on  Bode's  famous  gain-phase  integral  to  derive  the 
shape  of  a  nominal  loop  transfer  function.  The  method,  however,  needs  a  rational  function 
approximation  to  obtain  an  analytical  expression  for  the  loop  transfer  function,  and  straight  line 
approximations  for  the  nonlinear  QFT  bounds.  Ballance  and  Gawthrop  [1]  simplified  the 
iteration  process  of  this  technique,  using  the  QFT  Control  Design  Matlab  toolbox  [3]. 
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Thompson  and  Nwokah  [38]  proposed  a  method  based  on  nonlinear  programming  techniques 
wherein  the  templates  of  the  uncertain  plant  are  approximated  by  over  bounding  rectangles.  Such 
a  template  approximation  leads  to  over  bounding  in  the  constraints  derived  for  the  optimization. 

Bryant  and  Halikias  [5]  addressed  the  problem  of  automatic  loop  shaping  using  linear 
programming  techniques  wherein  the  QFT  bounds  are  approximated  by  a  series  of  linear 
approximations.  However,  their  method  also  leads  to  conservatism  in  describing  nonlinear  QFT 
bounds  by  linear  inequalities. 

Chait  et  al.  [6]  proposed  a  method  based  on  convexification  of  the  nonconvex  closed  loop 
bounds.  The  QFT  design  problem  in  this  method  is  posed  in  terms  of  the  closed  loop 
complementary  sensitivity  function.  The  novelty  of  this  method  is  that  the  closed  loop 
nonconvex  bounds  are  transformed  into  linear  inequalities  without  any  conservatism,  and  then  a 
linear  program  is  solved.  However,  as  pointed  by  these  authors,  the  shortcoming  of  the  method  is 
that  it  involves  fixing  the  poles  of  the  closed  loop  transfer  function  a  priori. 

The  main  drawback  of  the  approaches  given  above  lies  in  attempting  to  solve  a  complicated 
nonlinear  optimization  problem  using  convex  or  linear  programming  techniques,  which  generally 
leads  to  conservative  designs. 

The  section  1,  based  on  evolutionary  algorithms  [10]  and  the  section  2,  based  on  genetic 
algorithms  [11]  (as  Chen  et  al.  [9]),  will  try  to  overcome  these  difficulties,  reformulating  the 
problem  as  one  of  parameter  optimization  of  a  fixed  order  controller.  Those  are  stochastic 
methods,  which  are  based  on  a  probabilistic  approach.  Both  methods  are  computationally 
demanding,  and  it  is  well  known  that  with  genetic  and  evolutionary  algorithms  one  may  obtain  a 
local  minimum  instead  of  the  global  minimum.  However,  knowing  what  one  wants  helps  to  setup 
the  problem  in  a  way  that  increases  the  chance  of  getting  a  better  solution.  Besides,  it  has  the 
advantages  that  no  initial  controller  is  needed  and  that  it  uses  the  precise  values  of  numerical 
QFT  bounds,  which  overcomes  the  problems  associated  with  the  approximation  of  the  bounds 
mentioned  in  the  above  techniques. 

Finally,  Nataraj  and  Tharewall  [33]  proposed  an  interval  analysis  algorithm  to  automate  the 
controller  synthesis  step  of  QFT.  It  is  based  on  deterministic  interval  global  optimization 
techniques  and  also  uses  the  precise  values  of  numerical  QFT  bounds.  However,  the  controllers 
obtained  with  this  method  suffer  from  overdesign  over  the  design  frequency  range. 


2.2.  STRUCTURE 
2.2.1.  Introduction 

This  work  introduces  a  methodology  to  design  automatically  QFT  controllers  for  SISO  (single 
input-single  output)  plants  with  model  uncertainty.  The  method  generalizes  other  automatic  loop¬ 
shaping  techniques,  avoiding  restrictive  assumptions  about  the  search  space  -continuity, 
existence  of  derivatives,  convexity,  pre-specified  types  of  bounds,  and  other  matters-.  This 
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methodology  applies  the  Evolutionary  Algorithms  and  the  Genetic  Algorithms  to  search  the  QFT 
controller  that  fulfils  the  control  specifications  for  the  whole  set  of  models  with  uncertainty.  Both 
strategies  have  been  validated  with  several  benchmarks. 


2.2.2.  Methodology  description 

This  automatic  loop  shaping  methodology  has  been  applied  to  the  two-degree-of-freedom  control 
system  shown  in  Fig.l.  The  plant  under  consideration  P{s)  is  member  of  a  family  p,  exhibiting  a 
parametric  and  a  bounded  non-parametric  uncertainty  with  the  following  structure: 

p={p(a)[l  +  Aj:aG  a,A„G  A}  (1) 

where  \A^{jO)}<  m{o))  (2) 


D2(s) 


R(s) 


B 

Di(s) 


Y(s) 

-► 


N(s) 


Fig.  1.  Two-degree-of-freedom  feedback  structure 
with  external  disturbances  and  sensor  noise  signals. 

The  objective  is  to  design  the  pre-filter  F{s)  and  the  compensator  G{s)  so  that  specified  robust 
performance  is  achieved  despite  plant  uncertainty.  The  pre-filter  can  be  easily  designed  following 
the  methodology  in  Houpis,  Rasmussen  and  Garcia-Sanz  [30],  and  will  not  be  considered  here.  In 
this  work,  performance  specifications  are  not  limited  to  any  specific  type,  such  as  robust  stability, 
tracking,  disturbance  rejection,  etc.  In  fact,  any  type  is  allowed,  as  long  as  it  can  be  numerically 
computed. 

The  Fig. 2  shows  the  automatic  loop-shaping  process  proposed  in  this  work.  It  is  comprised  of  the 
following  steps: 

1.  -  Firstly  the  designer  has  to  define  the  system  plant  and  the  performance  specifications,  being 
then  translated  into  de  frequency  domain  and  obtaining  the  plant  templates  and  the  bounds 
respectively. 


On  the  other  hand,  the  user  has  to  define  the  search  algorithm  parameters.  Some  of  these 
parameters  are  the  next  ones:  -  Number  of  individuals  that  compose  the  population 

Number  of  generations 
Crossover  and  mutation  probability 
Fitness  function  weights  . . . 
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In  addition,  the  designer  has  to  define  the  structure  of  the  initial  controller  to  be  found.  Here  the 
designer  will  have  to  define  the  number  of  real  and/or  complex  poles  and  zeros,  as  well  as  its 
range  of  variation. 

2.  -  Select  the  strategy  to  search  the  optimal  controller.  In  this  work  the  Evolutionary  Algorithm 
and  the  Genetic  Algorithm  are  applied,  but  any  other  strategy  may  be  selected  in  this  step. 

3.  -  Once  the  search  algorithm  has  provided  the  best  solution,  it  is  saved  for  further  tasks. 

4.  -  Afterwards,  the  controller  obtained  with  the  search  algorithm  is  reduced  to  a  k-th  order 
model,  being  k  the  order  of  the  previous  controller  minus  1 .  This  is  performed  by  means  of  the 
optimal  Hankel  norm  approximation  [37].  The  resulting  model  is  saved  as  it  was  done  with  the 
previous  controller.  Then,  the  process  goes  back  to  the  step  number  2  and  calculates  the  next 
controller  with  a  structure  similar  to  the  one  resulted  from  applying  the  Hankel  norm 
approximation.  This  is  repeated  until  the  order  of  the  controller  is  equal  to  1. 

5.  -  Finally,  with  all  the  saved  results  some  tasks  may  be  performed  in  order  to  compare  them. 


Fig.  2.  Methodology  diagram 
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3.  AUTOMATIC  LOOP-SHAPING  OF  QFT  ROBUST  CONTROLLERS 
VIA  EVOLUTIONARY  ALGORITHMS 

3.1.  INTRODUCTION 

In  this  section  the  methodology  a)  applies  Evolutionary  Algorithms  to  search  the  QFT  controller 
that  fulfils  the  control  specifications  for  the  whole  set  of  models  with  uncertainty,  and  b) 
validates  the  new  approach  with  a  benchmark. 

Finally,  the  Matlab  files  are  described  in  Appendix  A. 


3.2.  STRATEGY  1:  EVOLUTIONARY  ALGORITHMS 


3.2.1.  Problem  Statement 


In  this  strategy  a  general  formulation  for  the  controller  structure  is  used.  It  is  expressed  by  the 
following  transfer  function: 


n,..  f  \  /  2j 

n 


r 


G{s)=k^ 


z=l 


\Zi  j 


/=1 


Vi^i 


-+ 


\  /  2 

n 


r 


j=i 


kPj  j 


j=i 


2 

Pi 

Pj 

(3) 


where,  kc  is  the  gain,  Zi  is  a  zero  that  may  be  complex  {ricz,  number  of  complex  zeros)  or  real  {rin, 
number  of  real  zeros),  and  pj  is  a  pole  (real  or  complex)  with  nirp  the  number  of  real  poles  and 
nicp  the  number  of  complex  poles.  Note  that  the  amount  of  complex  zeros  or  poles  must  be  even, 
to  have  pairs  of  complex  conjugate  numbers  and  obtain  a  polynomial  with  real  coefficients.  The 
controller  may  have  some  poles  in  the  origin  and  the  designer  can  check  parameter  r  (usually  0,  1 
or  2)  to  set  them. 


3.2.2.  Evolutionary  Algorithms 

Evolutionary  Algorithms  (FA)  ([23],  [32]  and  [34])  are  divided  in  three  major  branches: 
Evolutionary  Programming  (EP),  Evolution  Strategies  (ES)  and  Genetic  Algorithms  (GA).  Their 
common  features  are  the  use  of  a  population  and  the  exchange  of  information  among  individuals. 
They  apply  a  progressive  evolution  of  the  population  to  reach  a  final  individual  which  represents 
the  solution  of  the  problem.  These  algorithms  have  been  widely  applied  in  engineering  and 
operation  research,  particularly  on  numerical  optimization  problems.  These  applications  usually 
manage  several  objectives  subject  to  multiple  constraints,  so  that: 
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Minimize  [f^{x),  f^{x)] 

subject  to  g^{x)<0,g2{x)<0,...,g^{x)<0  (4) 

,X  =  {x,} 

In  multi-objective  optimization  [32]  the  concept  of  Pareto-optimal  solution  is  of  great 
importance.  A  possible  solution  Xk  is  called  Pareto-optimal  if  and  only  if  it  satisfies  the  following 
two  conditions: 

Vie  {1,2,..,, V) 

2S{1,2 . N] 

In  other  words,  Xk  is  Pareto-optimal  if  it  has  no  worse  values  for  all  the  objectives  than  any  other 
possible  solution  Xq,  and  at  least  it  improves  one  of  the  objectives.  The  problem  may  have  a  big 
number  of  Pareto-optimal  solutions,  and  it  could  be  impossible  to  find  all  of  them. 

Search  problems  are  often  established  to  look  for  as  many  solutions  as  possible.  However,  this 
strategy  deals  with  the  search  of  only  one  Pareto-optimal  solution.  Due  to  the  difficulty  of 
simultaneous  optimization  of  several  objectives  with  some  constraints,  the  multi-objective 
problem  is  reduced  to  a  single  objective  formulation  expressed  by  a  fitness  evaluation  function 
J{x)  with  a  trade-off  among  objectives  y,(v)  and  constraints  gj{x). 

N  M 

j{x)=Y, w,  fi {x)+^rjg  . {x)  (6) 

i=l  j=l 

Election  of  weights  (w„  rj)  is  critical  to  reach  a  good  performance.  This  is  usually  the  main 
handicap  of  the  single  function  approach.  Nevertheless,  alternatives  to  a  unified  fitness  function 
have  not,  at  the  moment,  clear  advantage  or  they  are  much  more  difficult  to  implement.  One  of 
the  reasons  for  the  popularity  of  evolutionary  algorithms  in  multi-objective  optimization  is  their 
simplicity  and  capacity  to  handle  with  complex  fitness  functions  and  non-linear  problems. 
However,  although  the  price  one  has  to  pay  is  a  big  computing  effort  in  comparison  with  other 
techniques,  that  is  perfectly  assumed  thanks  to  the  speed  and  capabilities  of  current  computers. 


3.2.3.  Proposed  Algorithm 

To  study  the  behaviour  of  evolutionary  algorithms  to  solve  the  QFT  controller  design  problem,  a 
simple  algorithm  based  on  GA  and  ES  has  been  developed.  The  main  features  that  define  this 
algorithm  are  described  in  the  present  section. 


1.3.3. 1.  Coding  Scheme 

A  real  coding  for  the  parameters  of  the  individuals  in  the  population  has  been  selected.  This  is  a 
characteristic  shared  with  the  ES,  in  contrast  with  the  classical  GA,  which  use  binary  strings. 
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Individuals  collect  all  controller  parameters  in  a  vector. 

x  =  [k^,z^,...,z„^,p^,...,p„^\  (7) 

where  =  tirz  +  ««  and  Up  =  riyp  +  ricp. 

In  case  of  a  complex  zero  or  pole,  only  the  absolute  value  of  the  real  and  imaginary  part  is 
considered.  Real  coding  avoids  coding-decoding  functions,  and  it  is  a  more  natural  approach  to  a 
numerical  optimization  of  engineering  problems. 


1.33.2.  Crossover 

Two  different  crossover  operators  have  been  used,  switching  between  one  and  another  randomly. 
In  particular,  the  discrete  crossover  (exchanging  parameters  between  individuals)  and  the 
intermediate  crossover  (weighting  parameters).  A  good  reference  of  them  may  be  found  in 
Michalewicz  [16].  The  definition  of  both  crossover  operators  follows  this  paragraph. 

Let  two  individuals  vi  =  ...,  Xn^)  and  V2  =  {xi'^,  ...,  Xn'^),  the  crossover  offspring  (vi2  = 

cross{v\,  vi))  may  be  obtained  by: 

•  Discrete  crossover:  vn  =  {xi  \  ...,  v„'") ,  where  i  =  { 1,2)  with  the  same  probability. 

•  Intermediate  crossover:  V2i  =  (a  xi'’^  +  (1-  a)  xi'’^,  ...  ,a  xj^  -i-  (1-  a  )  xj'^)  where  a  is  a  random 
number  between  0  and  1. 


1.3. 3. 3.  Mutation 

This  is  the  second  evolutionary  operator  to  create  new  individuals.  For  simplicity  a  uniform 
distribution  (a  bit  uncommon  in  ES)  is  defined,  which  varies  one  of  the  parameters  of  the 
individual  or  with  some  probability  all  of  them.  Therefore,  the  mutation  in  parameter  x,  of  vector 
V  may  be  expressed  as: 

mut{x, )  =  X,.  +  a{a  -  0.5){x,.  -  x,.  ^ )  (8) 

where  a  is  a  real  number,  so  that  0  <a  <  1 . 


1.3. 3.4.  Selection 

One  of  the  more  complex  issues  in  evolutionary  algorithms  is  the  selection  schema  of 
individuals.  In  this  sense,  a  deterministic  and  elitist  approach  is  applied,  where  population  is 
ordered  according  to  a  static  fitness  evaluation.  In  this  way,  a  population  of  p  individuals,  which 
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generates  X  offspring,  is  redueed  to  the  initial  size  seleeting  the  p  individuals  with  the  minimum 
fitness  funetion  from  the  total  p  +  A.  individuals. 

1.3.3. 5.  Algorithm  Procedure 

The  work  of  the  algorithm  may  be  summarized  in  the  following  steps. 

Step  1:  Create  an  initial  random  population  where  parameters  are  seleeted  into  their 
possible  ranges. 

Step  2:  Seleet  a  number  of  individuals  for  erossover  ('y\icross)  from  the  eurrent  population  (\|/d-  If 
the  number  is  not  even,  one  of  the  individuals  is  added  again.  The  erossover  is  performed 
sueeessively  every  two  individuals  ereating  a  new  population  as  a  result  (^cross_new)-  A  pair  of 
individuals  produees  two  offspring.  Considering  groups  of  four  individuals,  the  two  best  ones  are 
seleeted  to  survive  in  ^\icross_new 

Step  3:  A  number  of  individuals  are  seleeted  from  the  eurrent  population  to  perform  the 
mutation.  For  every  individual  the  offspring  is  ineluded  in  the  resulting  population  (\\imut)  if  it  is 
better  than  the  original  one.  The  non-mutated  individuals  are  ineluded  into  \\ino_mut- 

Step  4:  A  final  population  is  built  as  the  union  of  ^cross_new,  ^no_mut  and  The  size  is  redueed 
to  the  initial  one  seleeting  the  best  individuals.  They  are  the  population  for  the  next  generation 
(V?+r)- 

Step  5:  If  the  maximum  number  of  generations  has  not  been  reaehed,  jump  to  step  2. 


3.2.4.  Fitness  Evaluation 

As  it  was  mentioned  previously,  individuals  are  eompared  aeeording  to  the  evaluation  of  a  fitness 
funetion,  where  an  individual  is  better  than  any  other  if  it  eauses  a  lower  value  in  the  fitness 
evaluation.  This  funetion  is  defined  joining  in  an  expression  objeetives  and  eonstraints.  The 
eoeffieient  eleetion  is  very  eritieal  for  the  sueeess  in  the  seareh  proeess.  Definition  of  objeetives 
and  eonstraints  to  design  a  QFT  eontroller  is  a  subjeet  with  some  degrees  of  freedom.  It  has  been 
solved  in  the  following  manner. 


1.3.4. 1.  Objective  Functions 


The  set  of  objeetive  funetions  must  express  the  eontroller  quality  and  the  desired  behaviour  of  it 
on  the  NC.  To  simplify  the  mathematieal  expressions  in  the  rest  of  this  strategy,  the  following 
notation  is  introdueed: 


d{condition) 


1  ijf  condition  true 
0  ijf  condition  false 


(9) 
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The  set  of  work  frequencies  is  7,  and  sub-index  i  always  refers  to  a  frequency  of  this  set,  /e  7  . 
Also,  the  phase  of  Lo(s)  for  an  individual  v  is  expressed  as: 

^,,=Z4(x,ty,)  (10) 

To  improve  quality  a  minimization  of  gain  kc  is  required.  So,  the  first  objective  function /i(x)  can 
be  defined  as: 

fM=kl  (11) 

Two  other  objective  functions  are  included  to  manage  the  controller  behaviour  on  the  NC.  The 
function  f2{x)  express  an  heuristic  objective,  looking  for  controllers  with  a  Lo(s)  plot  as  near  as 
possible  to  the  UHFB  {Universal  High  Frequency  Boundary)  bound,  if  frequency  is  greater  than 
a  reference  COuhfb  (ty,  >  (^uhfb)-  This  reference  is  selected  among  the  set  of  frequencies  of  the 
QFT  analysis  where  bounds  start  to  be  in  convergence  with  the  UHFB.  Also,  the  phase  distance 
is  only  considered  if  it  is  bigger  than  the  reference  ^re/,  i.e.,  {\0uhfb  ~^xj  \  >  ^ref)^  where  (j) ref  is 
defined  by  the  designer.  See  the  next  Fig.3  to  clarify  the  content  of  this  paragraph. 


Fig.3.  Clarification  for  the  fiix)  objective  function. 


fl  (•^)  “  ^  ^  ^UHFB  )^{\^UHFB  ^x.i  \  ^  ^ ref  ^x,i 

i€l 


(12) 
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There  are  situations  where  it  is  very  useful  to  inelude  as  an  objeetive  the  seareh  of  one  eontroller 
aehieving  a  deseending  phase  for  Lo(s).  This  ean  be  set  by: 

(13) 


These  last  two  funetions  Ifiix),  fsix)]  might  not  be  eonvenient  in  some  situations.  Designer 
should  eheek  the  bounds  on  the  NC  to  deeide  whether  they  are  suitable.  They  are  ineluded  here 
beeause  of  the  usefulness  to  treat  some  speeifie  problems. 


1. 3.4.2.  Constraint  Functions 

Constraints  are  expressed  by  funetions  g ^  (v).  They  are  useful  to  know  whether  an  individual 

meets  a  eondition,  and  if  not,  to  measure  the  degree  of  mismateh.  When  an  individual  meets  all 
restrietions  it  is  ealled  a  feasible  individual.  Otherwise  it  is  an  unfeasible  solution  of  the  problem. 
In  the  next  paragraphs  the  funetions  that  have  been  eonsidered  for  this  study  are  introdueed. 
When  a  eontroller  meets  the  QFT  bounds  (over  upper  bounds  qu,  or  under  lower  bounds  qi,  see 
Fig. 4),  one  of  the  following  eonditions  is  true: 

,  s  I  \  IG  I  (14) 

Lo{x,O)J<qi[0^j,O)J 


Fig.4.  Clarification  for  Equation  (14). 

To  eheek  this  eondition,  and  quantify  the  level  of  bound  violation  if  it  is  neeessary,  the  following 
eonstraint  funetions  are  taken: 


V(x)='ZvUx.‘»,} 

iel 

and, 

g,{x)=S{v{x)>0){KM^)  +  '^) 


(15) 

(16) 
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where,  Ky=\lr^  with  the  weight  of  the  gi(x)  eonstraint  function,  and  where  the  sub¬ 
functions  V^{x,C0i)  explain  the  degree  of  mismatch  of  every  bound  at  frequency  <»,.  When  an 

individual  meets  their  correspondent  bounds,  the  value  is  0.  Otherwise  the  value  will  be  the 
minimum  between  the  distance  to  the  upper  and  to  the  lower  bound.  In  a  more  formal  way  it  can 
be  written  as: 

y„  {x,  CO. )  =  max{log  .  log  L„  {x,  CO, ),  o} 

V,  {x,  CO, )  =  maxjlog  Lq  {x,CO,)-\ogq,(e^,,CO,),0} 

and  finally, 

V  ix  CO  )  =  \  ) 

[max{y„  {x,  CO, ),  V,  {x,  O), )}  iff  q^  (o^,  ,0),)<  q,  (o^, ,  CO,  ) 

To  check  stability,  roots  of  the  characteristic  equation  [l-i-  Lo(v)]  are  computed.  The  real  part  of 
the  root  (J,)  is  the  distance  to  the  complex  axis.  Violation  of  the  stability  constraint  appears  if  the 
root  is  in  the  RHP  {Right-Half  Plane).  Therefore,  a  second  constraint  function  is  defined: 

4=X<J(</,  >oK 

and, 

g^{x)  =  d{dT  >0){K^,  d^+\) 

where,  K^,=\lr.^  with  the  weight  of  the  g2{x)  constraint  function. 

Finally,  sometimes  controllers  that  achieve  a  descending  module  plot  for  ILo(v)l  are  strongly 
preferred.  Every  controller  without  this  feature  is  considered  unfeasible.  This  avoids  controllers 
with  resonant  frequencies.  The  constraint  function  can  be  defined  as: 

g3 {x)  =  (v )|  always  descending]  (21) 

As  it  may  be  observed,  the  optimization  procedure  should  reach  feasible  solutions,  that 
means gi(v)  =  .?3(-^)  =  0,  trying  to  minimize  the  objective  functions /i(x), /2(v) 

and  /j  (x) .  These  functions  are  non-linear  and  present  a  rather  complex  behaviour. 


(19) 

(20) 


(17) 

(18) 


3.2.5.  Application  and  Results 

A  QFT  benchmark  example  is  proposed  to  test  the  quality  of  this  methodology.  This  benchmark  is 
the  same  as  the  one  selected  by  Chen  et  al.  [8],  Chait  et  al.  [6],  and  Borghesani  et  al.  [3]  in 
previous  works.  The  considered  plant  has  a  parameter  uncertainty  given  by: 
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P  = 


I  ke  [l,10],aG 


(22) 


The  nominal  plant  expressed  so  that  the  parameters  are  a=\  and  ^=1.  The  specifications 

for  robust  stability  and  reference  tracking  are  respectively: 


P{s)G{s) 
l  +  P(y)G(y) 


<1.2, 


VPg  P,ty > 0 


for  the  stability,  and. 


for  reference  tracking,  where. 


,_PA)Gi^ 

^  +  P(5)G(y) 


0.6854  (y  + 30) 
y'+4y  +  19.752 


and 


120 

+175^+825  +  120 


for  >  10  rad/s. 


Solutions  given  by  Chen  et  al.  [8]  and  Borghesani  et  al.  [3]  are: 


^Chen  (‘^) 


_ 1.7945  +  6.203 

0.0000045 '  +0.002745  +  0.826 


(23) 


(24) 


(25) 


(26) 


(27) 
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'  Borghesani 


is) 


3.0787xl0S" +3.5365x10^5+3.8529x10^ 
5^1.5288x10^5" +1.0636x10'’5 +4.281x10" 


(28) 


As  can  be  seen,  Chen  et  al.  [8]  provided  a  second  order  controller,  whereas  Borghesani  et  al.  [3] 
provided  a  third  order  one.  Their  characteristics  are  summarized  in  the  next  table: 


Lo(s)  (Chen) 

Lo(s)  (Borghesani) 

Phase  margin 

59.15° 

82° 

Gain  margin 

49.92 

52.66  dB 

HFG 

4.4852x10^ 

3.079x10^ 

Table  1:  Chen  and  Borghesani  controllers 


Note:  the  HFGs  (High  Frequency  Gain)  have  been  computed  according  to  the  formula  in  [6] : 

HFG=\ims^  G{s)  (29) 

where  r  is  the  excess  of  poles  of  Gf5). 

The  high-frequency  gain  (HFG)  of  4.4852x10^  in  Chen  is  smaller  than  that  of  Borghesani  of 
3.079X  10^.  This  means  that  the  Chen's  controller  needs  less  control  effort  and  is  less  sensitive  to 
high  frequency  noise.  In  addition,  while  the  gain  margins  are  quite  similar,  the  phase  margin  given 
by  Chen  (59.15°)  is  quite  smaller  than  the  one  given  by  Borghesani  (82°). 

Now,  the  methodology  developed  in  this  work  is  applied  to  the  problem  in  order  to  find 
controllers  and  to  compare  the  results.  To  begin  with,  the  same  controller  structure  of 
Borghesani 's  is  imposed:  a  gain,  two  real  zeros,  a  real  pole  and  a  pair  of  conjugate  complex  poles. 
Individual  is  expressed  by  the  search  vector  [kc,  zi,  Z2,  pi.reai,  Pumag,  pi]-  Wide  ranges  are  defined 
for  the  parameters:  e  [0.01,1500],  Zj  j  e  [0.01,400],  e  [0.01,1500],  e  [0.01,1500] 

and  P2  ^  [0.01,1500].  Throughout  the  execution,  the  controller  structure  changes  applying  the 
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optimal  Hankel  norm  approximation  [37]  to  obtain  a  model  reduction.  It  is  important  to 
emphasize  that  no  initial  controller  to  start  the  search  is  needed.  Only  an  initial  structure  and 
ranges  of  variation  for  the  parameters  are  needed.  That  is  one  of  the  strongest  features  of  the 
evolutionary  algorithms. 

The  complete  set  of  frequencies  to  loop  shape  the  controller  is  <2J=  {0.1,  0.5,  1,  2,  5,  10,  15,  40, 


10000)  (rad/s).  Crossover  probability  is  0.25,  mutation  probability  of  a  parameter  is  0.08  and 
a  =0.25 .  The  number  of  generations  is  300  and  the  population  size  is  80  individuals  in  each 
iteration.  Due  to  the  sensitivity  of  the  evolutionary  strategies  to  the  fitness  function  coefficients, 
some  of  them  have  been  considered  fixed  and  others  have  been  modified  for  different  simulations. 
Coefficients  and  K^t  are  0.25  (=l/r;_2),  =100 rad/s  and  Taking  this  into 

account,  the  methodology  achieves  the  following  solutions: 

Simulation  1: 


The  coefficients  for  the  fitness  evaluation  are  aj  =0.001,  aj  =0.001,  =  1,  <-,=[4,16], 

=  [4,16]  and  =  [4,16]. 

The  unique  feasible  solution  is  a  controller  with  the  same  structure  as  the  one  given  by  Borghesani 
et  al.  [3]: 


0.03009/ +3.935y+8.706 


(30) 


5.272x10"'/ +2.03x10"'/ +0.02685y+l 


It  is  quite  similar  to  the  one  proposed  by  Borghesani  et  al.  [3].  The  Nichols  plot  of  this  controller 
may  be  observed  in  Figure  7. 


60 


-360  -315  -270  -225  -180  -135  -90  -45  0 

Open-Loop  Phase  (deg) 


Fig.7.  Nichols  plot  for  Lo(s)  with  Gi(s) 
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According  to  this,  the  controllers  meet  the  problem  requirements.  The  loop  responses  in 
frequencies  corresponding  to  single-valued  bounds  are  above  their  respeetive  bound.  Upper 
frequeneies  are  out  of  the  UHFB  with  a  tangent  response. 

The  next  table  shows  the  some  performance  properties  of  this  eontroller: 


Lo(s)  (Gi) 

Phase  margin 

79.69° 

Gain  margin 

48.949  dB 

HFG 

5.7076x10^ 

Table  2:  Gi(s)  controller  performance  properties 


Simulation  2: 

The  eoeffieients  for  the  fitness  evaluation  are  aj  =0.001,  aj  =0.001,  a^=2,  rj=[7,25], 
Tj  =  [4,16]  and  =  [4,16]. 

The  eontrollers  obtained  in  this  simulation  make  Lo(s)  fulfil  the  bounds  at  every  frequency  defined 
for  the  bounds,  but  at  frequeneies  between  2  and  5  rad/s  the  Lo(s)  erosses  the  UHFB,  what  is  not 
desirable.  This  does  not  assure  the  elosed  loop  stability  at  those  frequencies  for  the  plant  with  its 
uncertainty.  The  solutions  are  interesting  even  though  they  are  not  feasible. 

The  controllers  are: 


_ 0.0105b^+3.092y  +  14.22 _ 

7.062  X  lO"""  y  V 1 .075  x  10"'  y "  +  0.005737  y  + 1 


(31) 
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3.0485  +  14.22 

2.049xl0“S" +0.0024135  +  1 


(32) 


Fig.9.  Nichols  plot  for  Lo(s)  with  Gsis) 


Fig. 8  and  Fig.9  show  that  the  same  loop-shape  is  obtained  with  both  controllers,  where  Gsis)  is 
the  reduced  model  of  Gsis)  after  applying  the  optimal  Hankel  norm  approximation  [32]. 


Uis)  iGs) 

Lois)  iGs) 

Phase  margin 

55.3° 

55.4° 

Gain  margin 

51.4  dB 

51.7  dB 

HFG 

1.4875x10^ 

1.4874x10^ 

Table  3:  G2(s)  and  GsCs)  controller’s  performance  properties 


Simulation  3: 

The  coefficients  for  the  fitness  evaluation  are  aj  =0.001,  aj  =0.001,  a^  =  2,  rj=[6,20], 
Tj  =  [4,16]  and  rj  =  [4,16]. 

The  controller  obtained  in  this  simulation  is: 

^  .  0.011515" +3.4545+7.522 

= -  (33) 

2.403x10“%" +9.379x10"®5" +0.012935 +1 

As  in  the  previous  simulation,  here  the  closed  loop  stability  is  not  assured  at  every  frequency  for 
the  whole  plant  uncertainty  range.  Specifically,  the  stability  is  not  assured  between  120  and  170 
rad/s  for  the  worst  case  within  the  uncertainty.  Apart  from  that,  all  the  bounds  are  fulfilled  at  the 
corresponding  frequencies. 
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Lo(s)  (G4) 

Phase  margin 

72.7° 

Gain  margin 

48.9  dB 

HFG 

4.79x10^ 

Table  4:  64(8)  controller  performance  properties 


Simulation  4: 

The  coefficients  for  the  fitness  evaluation  are  aj  =0.001,  6(2  =0.001,  a^=2,  rj=[5,18], 
Tj  =  [4,15]  and  Tj  =  [4,16]. 


The  feasible  solution  obtained  in  this  simulation  is: 

^  _ 0.02007y"+6.004y  + 8.499 _ 

“  9.772xl0“‘’y'+5.556xl0“S"+0.01173y  +  l 


(34) 


Fig.  11  shows  that  with  this  controller  all  the  performance  specification  bounds  are  fulfilled,  and 
that  the  closed  loop  stability  is  assured  at  every  frequency. 
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Lo(s)  (Gs) 

Phase  margin 

83.3° 

Gain  margin 

Al.l  dB 

HFG 

2.054x10^ 

Table  5:  Gs(s)  controller  performance  properties 


(Gchen) 

(^Borshesani) 

Lo(s)  (Gi) 

L„(s)  (G2) 

L„(s)  (Gs) 

Lo(s)  (G4) 

Lo(s)  (Gs) 

Phase 

margin 

59.15° 

0 

<N 

CO 

79.69° 

55.3° 

55.4° 

72.7° 

83.3° 

Gain 

margin 

49.92 

52.66  dB 

48.949  dB 

51.4  dB 

51.7  dB 

48.9  dB 

Al.l  dB 

HFG 

4.4852x10^ 

3.079x10'’ 

5.7076x10^ 

1.487x10'’ 

1.487x10'’ 

4.79X10^ 

2.054X10'’ 

Table  6:  All  the  controllers 


The  feasible  solutions  found  with  this  strategy  are  in  the  shadowed  columns.  Comparing  these 
feasible  solutions  (G;  and  Gj)  with  the  Chen's  and  Borghesani's,  it  is  clear  that  G/  and  Gj  are  quite 
similar  to  Gaen  and  GBorghesani  ■ 

The  G;  controller  is  quite  similar  to  Gsorghesani  in  terms  of  gain  and  phase  margin,  but  Gi  has  a 
smaller  HFG.  This  means  that  G;  needs  less  control  effort  and  is  less  sensitive  to  high  frequency 
noise.  Comparing  G;  with  Gchen  we  see  that  while  the  gain  margins  are  similar  and  the  HFG  are  of 
the  same  order,  the  phase  margin  is  lower  in  the  case  of  Gchen- 

The  Gj  controller  is  almost  equal  to  Gsorghesani  because  the  values  vary  a  little.  Comparing  it  with 
Gchen,  their  gain  margins  are  similar,  but  the  phase  margin,  as  well  as  the  HFG,  is  greater  in  the 
case  of  Gj  . 
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4.  AUTOMATIC  LOOP-SHAPING  OF  QFT  ROBUST  CONTROLLERS 
VIA  GENETIC  ALGORITHMS 

4.1.  INTRODUCTION 

In  this  section  the  methodology  a)  applies  Genetic  Algorithms  to  search  the  QFT  controller  that 
fulfils  the  control  specifications  for  the  whole  set  of  models  with  uncertainty,  and  b)  validates  the 
new  approach  by  using  several  benchmarks. 

The  corresponding  Matlab,  C-f-i-  and  Text  files  are  described  in  Appendix  B. 


4.2.  STRATEGY  2:  GENETIC  ALGORITHMS 
4.2.1.  Problem  Statement 

In  this  strategy  the  controller  structure  is  expressed  by  the  following  transfer  function: 

Y[Uio+Zi) 

G{x,jO))=k^^ -  (35) 

_p 

W{j(0+Ph) 

h=\ 

where,  x=[k(;  \  represents  the  parameters  in  the  search  space,  being  the 

gain,  Zi  the  real  zeros  dxxA  the  real  poles  (/i=l,...,np). 


4.2.2.  Genetic  Algorithms 

GAs  (Genetic  Algorithms)  are  search  procedures  based  on  the  mechanics  of  natural  selection  and 
natural  genetics.  A  population  of  candidate  solutions  (points  in  the  search  space)  is  maintained  at 
every  stage  or  every  generation,  using  GAs  nomenclature.  For  the  next  generation,  a  new  set  of 
artificial  creatures  is  produced  guided  by  the  survival  of  the  fittest  principle.  The  idea  is  to 
improve  performance  by  sampling  areas  of  the  search  space  that  are  more  likely  to  lead  to  good 
solutions.  Some  features  of  GAs  are  [23]: 

i.  The  parameter  set  x  is  coded  as  a  finite-length  string  over  some  finite  alphabet  F. 

ii.  GAs  search  from  a  population  of  such  strings,  reducing  the  probability  of  finding  local  peaks. 
Furthermore,  no  initial  point  is  used  to  start  from. 

iii.  GAs  use  only  the  information  provided  by  an  objective  (or  fitness)  function.  No  other 
auxiliary  information  is  required. 
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iv.  GAs  use  probabilistic  transition  rules  to  guide  the  search. 

A  simple  GA  is  composed  of  three  operators.  First,  reproduction,  by  which  individual  strings  are 
copied  to  the  next  generation.  Strings  with  a  higher  fitness  value  have  a  higher  probability  of 
contributing  one  or  more  offspring.  Then,  crossover,  by  which  some  newly  created  members  are 
mated  at  random.  From  two  mated  strings,  two  new  strings  are  created  by  swapping  some 
characters  between  them.  Finally,  mutation  is  the  random  alteration  (with  very  small  probability) 
of  the  value  of  a  character  in  the  string. 

The  mathematical  background  behind  GAs  is  the  concept  of  similarity  template  or  schema.  A 
similarity  template  is  a  string  from  the  extended  alphabet  F*  =  F  u  {*}  where  *  is  the  do  not  care 
symbol.  It  turns  out  that  the  number  of  such  templates  implicitly  processed  in  each  generation  is 
about  ,  where  n  is  the  number  of  members  in  the  population  (see  [23]  and  [25]).  This  implicit 
parallelism  gives  the  GA  much  of  its  power. 


4.2.3.  Proposed  Algorithm 

The  software  developed  for  the  GA-based  automatic  loop-shaping  is  a  C-i-i-  adaptation  of  the 
Simple  Genetic  Algorithm  in  [23].  The  design  decisions  made  in  the  program  are  the  following: 

i.  The  inverse  of  the  cost  function  (37)  is  used  -  GAs  are  usually  designed  to  maximise,  rather 
than  minimise. 

ii.  Fitness  scaling,  as  proposed  by  Goldberg  [23],  is  applied  to  obtain  the  actual  objective 
function. 


hi.  The  selected  alphabet  is  F  =  [0,  1],  that  is  to  say,  binary  strings  code  the  parameters  to 
optimise. 

iv.  The  i-th  parameter  in  the  search  space  is  coded  with  U  bits. 


V.  The  final  user  is  requested  to  provide  range  and  precision  of  the  parameters  by  entering  values 
for  Xi^min,  Xi^max  “  minimum  and  maximum  allowed  values  for  the  parameter  x,  -  and  for  the 
precision  7i„  defined  as. 


X  —  X 

7r.= - ; - ;  z  =  l,...,r;  r=l  +  n  +n 

2  '  -I 


(36) 


where  r  is  the  number  of  parameters,  a  value  the  user  can  also  modify.  Equation  (36)  is  used  by 
the  program  to  solve  for 

vi.  The  multi-parameter  coding  used  in  the  program  is  the  result  of  concatenating  r  single¬ 
parameter  codings.  Every  candidate  string  is  decoded  to  r  unsigned  integer  values  in  the  ranges 
[0,  2  ^‘■^],  and  then  these  values  are  mapped  linearly  to  the  ranges  [x,>m, 
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vii.  The  initial  population  is  randomly  generated. 

viii.  The  reproduction  operator  is  based  on  the  stochastic  remainder  selection  without 
replacement  (see  [4]  and  [23]),  where  the  expected  number  of  offspring  for  each  string  is 
calculated,  the  integer  part  of  these  values  directly  allocated  and  the  fractional  parts  treated  as 
probabilities  to  complete  the  population. 

ix.  Finally,  the  user  must  provide  the  following  additional  parameters  to  run  the  program: 

crossover  probability,  mutation  probability,  number  n  of  strings  in  every  population,  and  the  cost 
parameters  <», ,  ,  aj ,  aj  >  and  aj . 


4.2.4.  Fitness  Evaluation 

Following  the  cost  function  proposed  by  Thompson  and  Nwokah  [38]  and  Thompson  [39],  a  new 
reformulation  is  introduced  such  that. 


j[x)=aQk^  -l-aj  — — - a2A{o)i,0)2)+ a^V{x)  (37) 

,=i  PiZi 
®2 

where  (38) 

CO, 

and  where  V(x)  is  a  penalty  function,  to  be  discussed  below. 

This  new  cost  function  is  based  on  Horowitz’s  original  notion  of  optimum,  which  is  related  with 
the  well-known  cost  of  feedback  [26] .  The  area  term  allows  tighter  control  of  the  gain  over  any 
frequency  range  of  interest.  Indeed,  it  will  be  used  to  reduce  over-design  at  low  frequencies. 
Moreover,  additional  area  terms  for  other  frequency  ranges  could  be  added  when  necessary.  The 
V(x)  term  in  (37)  penalises  unfeasible  solutions,  i.e.  those  that  do  not  satisfy  the  stability  and 
performance  specifications. 

As  mentioned  above,  the  key  of  QFT  is  a  transformation  of  robust  stability  and  robust 
performance  specifications  and  plant  uncertainties  into  domains  in  the  complex  plane,  referred  to 
as  bounds,  where  a  nominal  loop  transmission  should  lie  within  [7].  A  generalised  bound  can  be 
represented  as  a  function  q  of  phase  and  frequency,  composed  of  upper  and  lower  parts,  qu  and  qi 
respectively,  such  that 


Lm  q^  (ZLq  [x,  jO). ),  0).  )<LmLQ  [x,  jO). ) 
Lm  L„  [x,  jO).  )<Lmq,  [ZL^  [x,  jO). ),  ) 


(39) 


where  Lo=GPo  and  Lm  denotes  logio  magnitude,  I  represents  the  discrete  set  of  frequencies  of 
interest  and  Po  is  the  selected  nominal  plant.  For  convenience,  it  is  denoted. 
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d{x,0).)=ZLQ{x,jCO.)  (40) 

[x,  CO. )  =  max  (Lm  [0[x,  CO. ),  CO. )  -  Lm  (x,  jco^  ),0  ) 

Vi  {x,  CO. )  =  max  (Lm  Lq  {x,  jco. )- Lm  q,  {0[x,  0)^ ),  CO^  ),0)  (41) 


represent  the  degree  of  violation  of  an  upper  or  lower  bound,  and  define 


min(y„  {x,co.),V,  {x,co.)), 
if  Lm  q^  (0(x,  CO. ),  CO.  )>Lmqi  (0(x,  CO^ ),  &»,. ), 
max(y„  (x,  CO.  ),V,  {x,  co^ )), 


otherwise. 


(42) 


Then,  the  penalty  function  is  defined  as  follows: 

vU)=X(V.0,®,)f  (43) 

iel 


which  completes  the  cost  function  in  (37). 

The  GA  program  simply  read  the  numerically  computed  bounds  from  a  text  file,  so  any  type  of 
bound  or  combinations  of  bounds  could  be  used.  For  convenience,  the  QFT  Toolbox  of 
MATLAB  [3]  is  used  to  generate  the  bounds.  Any  problem  that  can  be  stated  with  that  tool,  or 
with  any  other,  as  long  as  a  text  file  with  the  bounds  can  be  generated,  is  potentially  solvable 
with  this  methodology.  Also,  note  that  no  especial  requirements  are  imposed  on  the  plant 
templates. 

Typically,  CAD  tools  compute  bounds  only  for  a  discrete  set  of  phases  varying  in  a  certain  range, 
say  [-271,  0] .  In  order  to  evaluate  the  penalty  term  for  phase  values  where  qi  and  have  not  been 
computed,  a  linear  interpolation  was  used.  Of  course,  this  is  only  a  convenient,  simple-to- 
implement  solution.  Otherwise,  an  analytical  equation  for  the  bounds  could  be  directly 
programmed  when  available.  In  this  respect,  [7]  showed  that  there  exists  a  formal  map  from  the 
uncertain  plant  and  each  closed  loop  specification  to  the  bounds,  which  could  be  written  as 
programmable  quadratic  inequalities. 


4.2.5.  Application  and  Results 

Here  another  QFT  benchmark  example  is  proposed  to  test  the  quality  of  the  new  methodology. 
This  example  was  taken  by  Thompson  and  Nwokah  [38]  from  a  flight  control  problem  due  to 
Blakelock  [2].  The  plant  is: 
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1  +  - 

V  zj 


ft  0 

f 

2^  vM 

s 

1-I-  — 

1-1- 

1  pj 

V 

(Pn  0)1 

P{s)  = 


where  the  uncertain  parameters  are  given  by, 

[0.2,2],  zG  [0.5,75],  p& [l,10],  g  [5,6],  ^g  [0.8,0.9] 
The  nominal  plant  selected  for  the  loop-shaping  is 


1+- 

V 


1  +  - 

V  PoJ 


f 


^0=2,  Zo=0-5,  7*0  =10’  <^«o=0,  4=0-8 
Equations  (47),  (48),  (49)  and  (50)  exhibit  the  desired  closed-loop  specifications: 
P[s)G[s) 


1  +  P(v)g(v) 

for  the  stability,  and 


<6dfi,  VPg  P,  and  <»g  {0.01,0.05,0.1,0.2,1,5,10,50,100} 


F{s) 


p(s)a{s) 


1  +  P{s)g{s) 


for  reference  tracking,  where. 


Tjs)-- 


1  + 


^  V  ^ 

v0.35y 


( 

(  s\ 

1  +  — 

1  +  - 

1  0.5  j 

l  3j 

and 


4(^)  = 


(l  +  v)(l  +  v) 


7.T 

V  2y 


for  0)€  {0.01,0.05,0.1,0.2,1,5,10}. 


08-1-3027 
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(44) 


(45) 


(46) 


(47) 

(48) 

(49) 


(50) 
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Thompson  and  Nwokah  [38]  started  from  the  //^controller, 


7.709 


1  + 


V 


0.882 


1  + 


J\ 


2.22 


1  + 


J\ 


2.22 


1  +  - 

A  Vy 


1  +  - 

V  Vy 


1  + 


0.51 


1  + 


A 


0.615 


1  + 


A 


2.6 


y 


1  + 


28.36 


y 


1  + 


A 


1000 


(51) 


to  obtain,  after  a  non-linear  optimisation  procedure  and  limiting  each  controller  parameter  to  a 
40  percent  neighbourhood  of  its  initial  value, 

.  ^  37 1 .003  {s  +  0.630)(y  +  2.994)(y  +  2.994)(y  +  7.035)(y  +  7.035) 

™  ~  (y  +  0.7 14)(y  +  0.822)(y  +  3.640)(5  +  20.254)(y  +  7 14.29) 

As  required  by  their  optimisation  technique,  the  plant  templates  were  approximated  by  the 
smallest  possible  super  scribed  rectangle.  A  55%  reduction  in  the  high  frequency  gain  was 
obtained. 


Figure.  12  and  Table.7  show  the  loop-shaping  and  the  performance  properties  respectively  for  the 
Gjff  controller: 


Lo(s)  (Gtn) 

Phase  margin 

47.2° 

Gain  margin 

27.4  dB 

HFG 

3.71x10^ 

Table  7:  Gtn(s)  controller's  performance  properties 
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Note:  the  HFGs  (High  Frequency  Gain)  are  computed  according  to  the  formula  in  [6]: 

HFG  =  \\ms'^  G(v)  (53) 

where  r  is  the  excess  of  poles  of  G(s). 

Now,  the  GA  methodology  developed  in  this  work  is  applied  to  the  problem  in  order  to  compare 
and  validate  the  results.  To  begin  with,  the  same  controller  structure  of  Thompson  and  Nwokah  is 
imposed:  a  gain,  five  real  zeros  and  five  real  poles.  As  the  program  still  requires  maximum  and 
minimum  parameter  values,  wide  ranges  are  chosen:  k  g  [0.01,1000],  3  G  [0.01,5], 

£4  5  G  [0.01,10]  and  Pj  5  G  [0.01,1000].  With  such  wide  ranges,  the  GA  is  given  freedom  enough 
to  significantly  reduce  high  frequency  gain.  To  avoid  low  frequency  over-design,  the  area  term  in 
the  cost  function,  tuned  in  the  low  frequency  range,  is  used.  Cost  parameters  are  set  to  Gq  =0.001 , 
=0.00001,  a2=l,  area  of  frequency  range  [<»j  ,<»2]=[0.01,l],  and  ^3=!.  Throughout  the 
execution,  the  controller  structure  changes  applying  the  optimal  Hankel  norm  approximation  [37] 
to  obtain  a  model  reduction. 

The  complete  set  of  frequencies  to  loop  shape  the  controller  is  Q)=  (0.01  0.05  0.1  0.2  1  5  10  20 
30  40  50  100}  (rad/s).  Crossover  probability  is  0.1  and  the  mutation  probability  of  a  parameter  is 
0.01.  The  maximum  number  of  generations  is  500  and  the  population  size  is  100  individuals  in 
each  iteration.  Due  to  the  sensitivity  of  the  genetic  strategies  to  the  fitness  function  coefficients, 
some  of  them  have  been  considered  fixed  and  others  have  been  modified  for  different 
configurations. 

It  is  important  to  emphasize  that  no  initial  controller  to  start  the  search  is  needed.  Only  an  initial 
structure  and  ranges  of  variation  for  the  parameters  are  needed  (as  it  happened  with  the 
Evolutionary  Strategy  in  section  3). 

For  comparison  reasons  with  the  Thompson  and  Nwokah  [38]  controller,  the  rectangular 
approximation  for  the  templates  has  been  applied,  even  though  it  is  not  necessary  to  run  this 
methodology. 

As  it  was  mentioned  in  a  previous  paragraph,  different  configurations  are  chosen  varying  some 
parameters  while  others  are  kept  fixed.  The  results  obtained  are  the  following: 


Simulation  1: 


0.01 

ao  parameter 

0.0001 

ai  parameter 

1 

a2  parameter 

0.01 

Wini  to  measure  the 

"excess  controller  bandwidth" 

1 

wpin  to  measure  the 

"excess  controller  bandwidth" 

2 

a3  parameter 
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1440  gain  of  the  QFT  nominal  plant 

I  number  of  zeros  of  the  QFT  nominal  plant 

0.5  zero  of  the  nominal  plant ;  for  complex:  (2,3)  just  one  of  the  conjugates 

3  number  of  poles  of  the  QFT  nominal  plant  (without  counting  the  conjugates) 

0  pole  of  the  nominal  plant;  for  complex:  (2,3)  just  one  of  the  conjugates 

10  pole  of  the  nominal  plant;  for  complex:  (2,3)  just  one  of  the  conjugates 

(4.8,  3.6)  pole  of  the  nominal  plant;  for  complex:  (2,3)  just  one  of  the  conjugates 

I I  number  of  parameters 

0.01  minimum  value  for  the  parameter  1 
0.01  minimum  value  for  the  parameter  2 
0.01  minimum  value  for  the  parameter  3 
0.01  minimum  value  for  the  parameter  4 
0.01  minimum  value  for  the  parameter  5 
0.01  minimum  value  for  the  parameter  6 
0.01  minimum  value  for  the  parameter  7 
0.01  minimum  value  for  the  parameter  8 
0.01  minimum  value  for  the  parameter  9 
0.01  minimum  value  for  the  parameter  10 
0.01  minimum  value  for  the  parameter  1 1 
1000  maximum  value  for  the  parameter  1 
5  maximum  value  for  the  parameter  2 

5  maximum  value  for  the  parameter  3 

5  maximum  value  for  the  parameter  4 
10  maximum  value  for  the  parameter  5 
10  maximum  value  for  the  parameter  6 
1000  maximum  value  for  the  parameter  7 
1000  maximum  value  for  the  parameter  8 
1000  maximum  value  for  the  parameter  9 
1000  maximum  value  for  the  parameter  10 
1000  maximum  value  for  the  parameter  1 1 
.01  desired  precision  for  parameters 
200  number  of  individuals 

0. 1  crossover  probability 

0.01  mutation  probability 

1  show  continuously  (1)  or  just  at  the  end  (0) 

500  maximum  number  of  generations 
100000  minimum  desired  fitness 

0  number  of  controller  fixed  poles 

5  number  of  parameters  that  corresponds  to  controller  zeros 

The  controllers  obtained  with  this  configuration  are: 
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G,{s) 


126.71  l(y  +  3.31  l)(y  +  3.74)(y  +  4.053)(y  +  4.16)(y  +  5.81 1) 
{s  +  1.017)(y  +  1.063)(y  +  5.93)(y  +  23.19)(y  +  224.2) 


(54) 


for  the  original  structure  (similar  to  Gtn)- 


G^is) 


126.71 1  (y  ^  +  6.476y  +  10.63)(y "  +  8.672y  + 19.27) 
{s  +  224.2){s  +  23.19)(y  +  1.07)(y  + 1.01 1) 


(55) 


with  a  structure  simplified  once. 
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Lo(s)  (Gi) 

Lo(s)  (G2) 

Phase  margin 

52.9° 

52.9° 

Gain  margin 

19.2  dB 

19.2  dB 

HFG 

1.26711x10^ 

1.26711x10^ 

Table  8:  Gi(s)  and  G2(s)  controllers'  performance  properties 


Simulation  2: 


0.01  ao  parameter 

0.0001  ai  parameter 

1  a2  parameter 

0.01  Wini  to  measure  the  "excess  controller  bandwidth" 

1  Wpin  to  measure  the  "excess  controller  bandwidth" 

2  a3  parameter 

1440  gain  of  the  QFT  nominal  plant 
1  number  of  zeros  of  the  QFT  nominal  plant 

0.5  zero  of  the  nominal  plant ;  for  complex:  (2,3)  just  one  of  the  conjugates 

3  number  of  poles  of  the  QFT  nominal  plant  (without  counting  the  conjugates) 

0  pole  of  the  nominal  plant;  for  complex:  (2,3)  just  one  of  the  conjugates 

10  pole  of  the  nominal  plant;  for  complex:  (2,3)  just  one  of  the  conjugates 

(4.8,  3.6)  pole  of  the  nominal  plant;  for  complex:  (2,3)  just  one  of  the  conjugates 


1 1  number  of  parameters 
0.01  minimum  value  for  the  parameter  1 
0.01  minimum  value  for  the  parameter  2 
0.01  minimum  value  for  the  parameter  3 
0.01  minimum  value  for  the  parameter  4 
0.01  minimum  value  for  the  parameter  5 
0.01  minimum  value  for  the  parameter  6 
0.01  minimum  value  for  the  parameter  7 
0.01  minimum  value  for  the  parameter  8 
0.01  minimum  value  for  the  parameter  9 
0.01  minimum  value  for  the  parameter  10 
0.01  minimum  value  for  the  parameter  1 1 
1000  maximum  value  for  the  parameter  1 
5  maximum  value  for  the  parameter  2 
5  maximum  value  for  the  parameter  3 
5  maximum  value  for  the  parameter  4 
10  maximum  value  for  the  parameter  5 
10  maximum  value  for  the  parameter  6 
1000  maximum  value  for  the  parameter  7 
1000  maximum  value  for  the  parameter  8 
1000  maximum  value  for  the  parameter  9 
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1000  maximum  value  for  the  parameter  10 

1000  maximum  value  for  the  parameter  1 1 

.01  desired  precision  for  parameters 

300  number  of  individuals 

0.1  crossover  probability 

0.01  mutation  probability 

1  show  continuously  (1)  or  just  at  the  end  (0) 

500  maximum  number  of  generations 
100000  minimum  desired  fitness 

0  number  of  controller  fixed  poles 

5  number  of  parameters  that  corresponds  to  controller  zeros 


The  controllers  obtained  with  this  configuration  are: 


105.662  (y  +  4.063)(y  +  4.385)(y  +  4.385)(y  +  4.883)(y  +  9.092) 
{s  +  16.1l){s  + 1 1 .42)(5  +  159.4)(y  + 1 .566)(5  + 1 .322) 


for  the  original  structure  (similar  to  Gtn)- 


(56) 


G4(^) 


105.662  (/  +5.677y  + 10.05 )(/  +10.67y  + 36.37) 
{s  +  157.4)(5  +  28.08)(5  +  \.91A){s  +  1.217) 


with  a  structure  simplified  once. 


(57) 
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G,{s) 


63.0897  (y  +  4.96  l)(y  +  4.365)(y  +  3.525)(y  +  8.164) 
{s  +  0.4494)(5  +  4.328)(5  +  43.46)(5  +  71.38) 


(58) 


applying  the  GAs  with  the  previous  simplified  strueture. 


Lo(s)  (G3) 

Lo(s)  (G4) 

Lo(s)  (Gs) 

Phase  margin 

53.2° 

53.4° 

49.2° 

Gain  margin 

16.6  dB 

16.6  dB 

12.8  dB 

HFG 

1.05662x10^ 

1.05662x10^ 

0.6309x10^ 

Table  9:  GsCs),  64(8)  and  Gs(s)  controllers'  performance  properties 
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Simulation  3: 

0.01  ao  parameter 

0.0001  ai  parameter 

1  a2  parameter 

0.01  Wini  to  measure  the  "excess  controller  bandwidth" 

1  Wpin  to  measure  the  "excess  controller  bandwidth" 

3  a3  parameter 

1440  gain  of  the  QFT  nominal  plant 

I  number  of  zeros  of  the  QFT  nominal  plant 

0.5  zero  of  the  nominal  plant ;  for  complex:  (2,3)  just  one  of  the  conjugates 

3  number  of  poles  of  the  QFT  nominal  plant  (without  counting  the  conjugates) 

0  pole  of  the  nominal  plant;  for  complex:  (2,3)  just  one  of  the  conjugates 

10  pole  of  the  nominal  plant;  for  complex:  (2,3)  just  one  of  the  conjugates 

(4.8,  3.6)  pole  of  the  nominal  plant;  for  complex:  (2,3)  just  one  of  the  conjugates 

I I  number  of  parameters 

0.01  minimum  value  for  the  parameter  1 
0.01  minimum  value  for  the  parameter  2 
0.01  minimum  value  for  the  parameter  3 
0.01  minimum  value  for  the  parameter  4 
0.01  minimum  value  for  the  parameter  5 
0.01  minimum  value  for  the  parameter  6 
0.01  minimum  value  for  the  parameter  7 
0.01  minimum  value  for  the  parameter  8 
0.01  minimum  value  for  the  parameter  9 
0.01  minimum  value  for  the  parameter  10 
0.01  minimum  value  for  the  parameter  1 1 
1000  maximum  value  for  the  parameter  1 
5  maximum  value  for  the  parameter  2 
5  maximum  value  for  the  parameter  3 
5  maximum  value  for  the  parameter  4 

10  maximum  value  for  the  parameter  5 
10  maximum  value  for  the  parameter  6 
1000  maximum  value  for  the  parameter  7 
1000  maximum  value  for  the  parameter  8 
1000  maximum  value  for  the  parameter  9 
1000  maximum  value  for  the  parameter  10 
1000  maximum  value  for  the  parameter  1 1 
.01  desired  precision  for  parameters 
250  number  of  individuals 
0.1  crossover  probability 

0.01  mutation  probability 
1  show  continuously  (1)  or  just  at  the  end  (0) 

500  maximum  number  of  generations 
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100000  minimum  desired  fitness 
0  number  of  controller  fixed  poles 

5  number  of  parameters  that  corresponds  to  controller  zeros 


The  controllers  obtained  with  this  configuration  are: 

^  ^  _  282.968 {s  +  4.023)(y  +  4.463)(y  +  4.902)(y  +  5.78 l)(y  +  6.787) 

®  “  (y  +  0.3533)(y  +  2.985)(5  +  7.327 Xv  +  52.2l){s  +  249.7) 


(59) 


for  the  original  structure  (similar  to  Gtn)- 


G,{s) 


282.968  (y  ^  +  7.61y  + 15.07  )(y"  + 1 1. 15y  +  32.58) 
(y  +  249.6)(y  +  52.3  l)(y  +  3. 107)(y  +  0.3533) 


with  a  structure  simplified  once. 


(60) 
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^  r  48.777 (y  +  3.477)(y  +  3.818)(y  +  3.916)(y  + 4.356) 
“  (y  + 1 . 104)(5  +  0.7806)(5  +  30.95)(5  +  70.96) 


(61) 


applying  the  GAs  with  the  previous  simplified  structure. 


Lo(s)  (Ge) 

Lo(s)  (G7) 

Lo(s)  (Gs) 

Phase  margin 

64.4° 

64.4° 

47.8° 

Gain  margin 

20.5  dB 

20.5  dB 

12.1  dB 

HFG 

2.829x10^ 

2.829x10^ 

0.487x10^ 

Table  10:  G6(s),  67(8)  and  GgCs)  controllers'  performance  properties 
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Simulation  4: 

0.01  ao  parameter 

0.0001  ai  parameter 

1  a2  parameter 

0.01  Wini  to  measure  the  "excess  controller  bandwidth" 

1  Wpin  to  measure  the  "excess  controller  bandwidth" 

3  a3  parameter 

1440  gain  of  the  QFT  nominal  plant 

I  number  of  zeros  of  the  QFT  nominal  plant 

0.5  zero  of  the  nominal  plant ;  for  complex:  (2,3)  just  one  of  the  conjugates 

3  number  of  poles  of  the  QFT  nominal  plant  (without  counting  the  conjugates) 

0  pole  of  the  nominal  plant;  for  complex:  (2,3)  just  one  of  the  conjugates 

10  pole  of  the  nominal  plant;  for  complex:  (2,3)  just  one  of  the  conjugates 

(4.8,  3.6)  pole  of  the  nominal  plant;  for  complex:  (2,3)  just  one  of  the  conjugates 

I I  number  of  parameters 

0.01  minimum  value  for  the  parameter  1 
0.01  minimum  value  for  the  parameter  2 
0.01  minimum  value  for  the  parameter  3 
0.01  minimum  value  for  the  parameter  4 
0.01  minimum  value  for  the  parameter  5 
0.01  minimum  value  for  the  parameter  6 
0.01  minimum  value  for  the  parameter  7 
0.01  minimum  value  for  the  parameter  8 
0.01  minimum  value  for  the  parameter  9 
0.01  minimum  value  for  the  parameter  10 
0.01  minimum  value  for  the  parameter  1 1 
1000  maximum  value  for  the  parameter  1 
5  maximum  value  for  the  parameter  2 
5  maximum  value  for  the  parameter  3 
5  maximum  value  for  the  parameter  4 

10  maximum  value  for  the  parameter  5 
10  maximum  value  for  the  parameter  6 
1000  maximum  value  for  the  parameter  7 
1000  maximum  value  for  the  parameter  8 
1000  maximum  value  for  the  parameter  9 
1000  maximum  value  for  the  parameter  10 
1000  maximum  value  for  the  parameter  1 1 
.01  desired  precision  for  parameters 
300  number  of  individuals 
0.2  crossover  probability 

0.01  mutation  probability 
1  show  continuously  (1)  or  just  at  the  end  (0) 

500  maximum  number  of  generations 
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100000  minimum  desired  fitness 
0  number  of  controller  fixed  poles 

5  number  of  parameters  that  corresponds  to  controller  zeros 


The  controllers  obtained  with  this  configuration  are: 

^  ^  123.949 (y  + 2. l)(y  +  4.063)(y  +  4.6)(y  + 9.805) 

“  (y  +  0.879l)(y  +  1.528)(y  +  81.18)(5  +  92.l) 

applying  the  GAs  with  a  structure  simplified  once. 


Lo(s)  (G9) 

Phase  margin 

63.4° 

Gain  margin 

17.4  dB 

HFG 

1.239x10^ 

Table  11:  69(8)  controller’s  performance  properties 


Simulation  5: 


0.1  ao  parameter 

0.0001  ai  parameter 
1  a2  parameter 

0.01  Wini  to  measure  the  "excess  controller  bandwidth" 

1  wpin  to  measure  the  "excess  controller  bandwidth" 

3  aa  parameter 
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1440  gain  of  the  QFT  nominal  plant 

I  number  of  zeros  of  the  QFT  nominal  plant 

0.5  zero  of  the  nominal  plant ;  for  complex:  (2,3)  just  one  of  the  conjugates 

3  number  of  poles  of  the  QFT  nominal  plant  (without  counting  the  conjugates) 

0  pole  of  the  nominal  plant;  for  complex:  (2,3)  just  one  of  the  conjugates 

10  pole  of  the  nominal  plant;  for  complex:  (2,3)  just  one  of  the  conjugates 

(4.8,  3.6)  pole  of  the  nominal  plant;  for  complex:  (2,3)  just  one  of  the  conjugates 

I I  number  of  parameters 

0.01  minimum  value  for  the  parameter  1 
0.01  minimum  value  for  the  parameter  2 
0.01  minimum  value  for  the  parameter  3 
0.01  minimum  value  for  the  parameter  4 
0.01  minimum  value  for  the  parameter  5 
0.01  minimum  value  for  the  parameter  6 
0.01  minimum  value  for  the  parameter  7 
0.01  minimum  value  for  the  parameter  8 
0.01  minimum  value  for  the  parameter  9 
0.01  minimum  value  for  the  parameter  10 
0.01  minimum  value  for  the  parameter  1 1 
1000  maximum  value  for  the  parameter  1 
5  maximum  value  for  the  parameter  2 
5  maximum  value  for  the  parameter  3 
5  maximum  value  for  the  parameter  4 
10  maximum  value  for  the  parameter  5 
10  maximum  value  for  the  parameter  6 
1000  maximum  value  for  the  parameter  7 
1000  maximum  value  for  the  parameter  8 
1000  maximum  value  for  the  parameter  9 
1000  maximum  value  for  the  parameter  10 
1000  maximum  value  for  the  parameter  1 1 
.01  desired  precision  for  parameters 

250  number  of  individuals 

0.15  crossover  probability 

0.01  mutation  probability 

1  show  continuously  (1)  or  just  at  the  end  (0) 

500  maximum  number  of  generations 

100000  minimum  desired  fitness 

0  number  of  controller  fixed  poles 

5  number  of  parameters  that  corresponds  to  controller  zeros 

The  controllers  obtained  with  this  configuration  are: 
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gM 


62.7769  (y  +  2.705 )(y  +  4. 17)(y  +  4.707 )(y  +  5.0l)(y  +  8.69 1) 
{s  + 1.71  l)(y  +  0.5517)(y  +  7.639)(y  +  44.85)(y  +  69.92) 


(63) 


for  the  original  structure  (similar  to  Gtn)- 


gM 


62.7769(y  +  7.283)(y  +  2.39l)(y^  +7.845^  +  16.79) 
{s  +  70.62)(5  +  44.24)(y  + 1.657)(5  +  0.5528) 


(64) 


with  a  structure  simplified  once. 


41 


Advanced  Quantitative  Robust  Control  Engineering:  Ref.  BOARD  Award  #  FA8655-08-1-3027 

New  Solutions  for  Automatic  Loop-shaping  for  SISO  Ref.  Grant  #  083027 

and  MIMO  Systems.Part  1:  SISO  Systems.  Date:  1  .September.2009  Page  42 


Lo(s)  (Gio) 

Lo(s)  (Gn) 

Phase  margin 

52.4° 

52.3° 

Gain  margin 

13.2  dB 

13.2  dB 

HFG 

0.627x10^ 

0.627x10^ 

Table  12:  Gio(s)  and  Gii(s)  controllers'  performance  properties 


L„(s) 

L„(s) 

L„(s) 

L„(s) 

L„(s) 

L„(s) 

Lo(s) 

Lo(s) 

Lo(s) 

Lo(s) 

L„(s) 

Lo(s) 

(G™) 

(Gd 

(G2) 

(Gs) 

(G4) 

(Gs) 

(GJ 

(Gy) 

(Gs) 

(G,) 

(Gio) 

(Gn) 

Phase 

margin 

A12° 

52.9° 

52.9° 

53.2° 

53.4° 

49.2° 

64.4° 

64.4° 

47.8° 

63.4° 

52.4° 

52.3° 

Gain 

21 A 

19.2 

19.2 

16.6 

16.6 

12.8 

20.5 

20.5 

12.1 

17.4 

13.2 

13.2 

margin 

dB 

dB 

dB 

dB 

dB 

dB 

dB 

dB 

dB 

dB 

dB 

dB 

HFG 

371 

126.7 

126.7 

105.6 

105.6 

63.09 

282.9 

282.9 

48.7 

123.9 

62.7 

62.7 

Number 
of  poles 

10 

10 

8 

10 

8 

8 

10 

8 

8 

8 

10 

8 

+  zeros 

Table  13:  All  controllers'  performance  properties 


Table  13  compares  the  achieved  controllers  with  the  one  given  by  Thompson  and  Nwokah  [38]. 
As  we  can  see,  the  phase  margin  of  the  new  controllers  is  always  higher  (but  close  to  47.2°),  the 
gain  margin  is  always  lower  and  the  High  Frequency  Gain  is  in  most  cases  quite  lower  than  37 1 . 
In  addition,  the  new  controllers  have  a  simpler  (or  the  same)  structure  than  the  Thompson  and 
Nwokah  one.  With  these  results  it  is  possible  to  conclude  that  the  new  GA  methodology  applied 
in  this  work  is  good  enough  to  achieve  useful  controllers  as  a  first  approach. 
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5.  CONCLUSION 

As  it  was  concluded  with  the  Evolutionary  Algorithms,  the  results  obtained  with  the  Genetic 
Algorithm  show  that  a  good  feasible  solution  is  possible  to  achieve  with  the  new  automatic  loop¬ 
shaping  strategy.  At  the  same  time  it  is  important  to  accentuate  that  the  strategy  is  sensitive  to  the 
fitness  function  coefficients.  Depending  on  their  adjustment,  a  feasible  or  an  unfeasible  solution 
is  achieved.  In  that  point,  it  is  possible  to  conclude  that  intuitive  considerations  in  the  search  may 
be  very  useful,  and  the  adjustment  of  the  fitness  evaluation  function  is  of  great  importance  for  the 
success.  Defining  functions  less  sensitive  to  coefficient  value  should  be  an  interesting  future 
work  because  this  is,  in  fact,  one  of  the  main  limitations  of  the  Evolutionary  and  Genetic 
Algorithms. 

Although  not  always  the  designer  can  achieve  a  feasible  solution  to  the  problem,  it  can  represent 
a  first  approach  to  the  final  solution,  in  the  sense  that  the  obtained  controller  may  be  modified 
until  getting  a  feasible  one. 

To  sum  up,  this  methodology  based  on  Evolutionary  Algorithms  and  Genetic  Algorithms  can 
contribute  to  improve  QET  controller  design,  alleviating  much  of  the  manual  work  required  at 
present. 
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7.  APPENDIX  A:  MATLAB  FILES  FOR  THE  STRATEGY  I 
(EVOLUTIONARY  ALGORITHMS) 


%  AUTOLOOP SHAPING_ES  Program  to  synthesize  optimal  controllers  with  an 
%  initial  structure  which  can  be  modified  by  means  of 

%  optimal  Hankel  norm  approximation  If  any  modification 

%  occurs  the  program  synthesizes  the  optimal  controllerf or 

%  the  new  structure. 

o, 

0 

%  Authors:  Mario  Garcia-Sanz,  Carlos  Molins 
%  1 . September . 2 00 9 


% - Plant  templates  and  bounds - 

%  Before  executing  this  program  the  user  has  to  define  the  plant  and  the 
%  specifications  at  Tmpl_Bnds.m,  as  well  as  all  the  parameters  needed  for 
%  the  search  strategy  at  configuration . m 

% - Templates  and  bounds  creation - 

[bdb,  np,  dp]  =  Tmpl_Bnds ( )  ; 

% - Parameters  configuration - 

%  Configuration  of  all  the  parameters  needed  for  the  search  strategy, 
configuration; 

% - Automatic  loop-shaping - 

z  =  l; 

nr z  0=nr z  ; 
ncz  0=ncz  ; 
npr0=npr; 
npc0=npc; 
k=npr+2  *npc-l  ; 

while  k>=l  %  Loop  to  calculate  the  controllers  while  the  structure  varies. 
% - Evolutive  Algorithm - 


EvAl  ; 

%  The  results  obtained  in  each  iteration  are  stored  in  vectors 
%  (gain,  realzeros, complexzeros  .  .  . ) 

gain (z, : ) =kg; 

if  nrz==0 

realzeros(z, :)=0; 

else 

realzeros (z, 1 :nrz) =rz; 

end 
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if  ncz==0 

complexzeros ( z  ,  :  ) =0  ; 

else 

complexzeros (z, 1 :2*ncz) =cz; 

end 

if  npr==0 

realpoles  ( z ,  : ) =0 ; 

else 

realpoles ( z , 1 : npr ) =pr; 

end 

if  npc==0 

complexpoles ( z , : ) =0 ; 

else 

complexpoles ( z , 1 : 2  *npc) =pc; 

end 

f cost ( z , : ) =y ; 
z=z+l ; 


% - Model  reduction - 

o, 

0 

%  Now  the  obtained  controller  is  reduced  to  a  k-th  order  model  by  means 
%  of  the  optimal  Hankel  norm  approximation. 

[numC,  denC]  =  TG (kg, cz , rz , pc, pr ) ; 

sys=tf (numC, denC) ; 

sys=ss (sys) ; 

n=size (sys. A, 1) ; 

nu=size (sys, 2) ; 

[sysc,  cinfo] =ncfmr (sys, n) ;  %  or  obtain  coprime  factors  of  system 

sysch=hankelmr (cinf o . GL, k) ;  %  or  optimal  Hankel  norm  approximation 
syschm=minreal ( inv ( sysch ( : , nu+1  ...%and  obtain  k-th  order  model 
: end) ) *sysch (:,l:nu)); 

[ Z , P , K] =ZPKDATA ( syschm, ' V ' ) ;  %  Gain,  poles  and  zeros  of  the  reduced 

%  model. 

%  Storage  of  the  gain,  poles  and  zeros  of  the  reduced  model  in  the 
%  correspondent  vectors  (gain,  complexzeros,  realzeros, . . .) . 

n  r  z  =  0  ; 
ncz  =  0  ; 
npr=0 ; 
npc=0 ; 

realzeros (z, 1 :nrz) =0; 
complexzeros (z,l:2*ncz)=0; 
realpoles ( z , 1 : npr ) =0 ; 
complexpoles ( z , 1 : 2  *npc) =0 ; 

if  ~isempty(Z) 

11=1; 
ri=l  ; 
ci=l  ; 

while  ii<=length ( Z ) 

if  imag ( Z ( 11 ) ) ~=0 

complexzeros ( z , cl ) =abs (real ( Z ( 11 )  )  )  ; 
complexzeros ( z , ci+1 ) =abs ( imag ( Z ( 11 ) ) ) ; 
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ii=ii+2 ; 
ci=ci+2 ; 
ncz=ncz+l ; 

else 

realzeros (z, ri) =abs  (Z  (ii) ) ; 
ii=ii+l ; 
ri=ri+l ; 
nrz=nrz+l ; 


end 

end 

end 


if  ~isempty(P) 

j  j  =  i; 

r  j  =  l; 
c  j  =  l; 

while  j j<=length (P) 

if  imag (P ( j  j ) ) ~=0 

complexpoles { z , c j ) =abs { real (P { j  j  )  )  )  ; 
complexpoles { z , c j  +  1 ) =abs { imag (P { j  j )  )  )  ; 
j j=j j+2; 

c  j=c  j+2  ; 
npc=npc+l ; 

else 

realpoles ( z , r j ) =abs (P { j  j ) )  ; 

j j=j j+1; 

r j=r j+1 ; 
npr=npr+l ; 

end 

end 

end 

kg=K; 

rz=realzeros (z, l:nrz) ; 
cz=complexzeros (z,l:2*ncz); 
pr=realpoles ( z , 1 : npr ) ; 
pc=complexpoles ( z , 1 : 2 *npc)  ; 

Kg  =  K_cor2 (kg, cz , rz , pc,  pr)  ; 
gain  (z,  : ) =Kg; 

z  =  z  +  l ; 

k=npr+2  *npc-l  ; 

Range_z (:,  (nrz+2*ncz  +  l)  : end)  =  [ ] ; 

Range_p ( : ,  (npr+2  *npc+l )  : end)  =  [ ] ; 

if  k==0  %  This  loop  is  to  calculate  the  1st  order  controller  with  the 
%  search  strategy. 

EvAl  ; 

gain  (z,  : ) =kg; 

if  nrz==0 

realzeros(z, :)=0; 

else 

realzeros (z, 1 :nrz) =rz; 


end 

if  ncz==0 


complexzeros ( z ,  : ) =0  ; 

else 
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complexzeros (z,l:2*ncz)=cz; 

end 

if  npr==0 

realpoles  ( z ,  : ) =0 ; 

else 

realpoles ( z , 1 : npr ) =pr  ; 

end 

if  npc==0 

complexpoles ( z , : ) =0 ; 

else 

complexpoles ( z , 1 : 2  *npc) =pc; 

end 

f cost ( z , : ) =y ; 
z=z+l ; 

end 

end 


%  Once  the  controllers  have  been  obtained  some  optional  tasks  can  be 
%  performed 

clear  cp; 

[r,  t] =size (np) ; 
for  i=l : 1 : r 

cp  (1,  1, 1)  =tf  (np  (i,  :  )  ,  dp  (i,  :  )  )  ; 

end 


%  Storage  of  all  the  obtained  controllers 
controllers ; 

%  Menu  to  perform  several  tasks  with  the  obtained  controllers. 

wcheck=logspace (-4, 4,  100)  ;  save  wcheck;  %  frequencies  for  the  chksiso  in 
menuplot 

wloop=logspace (-2, 4, 200) ;  save  wloop;%  frequencies  for  the  Ipshape  in  menuplot 
menuplot ; 


%  BVALID  Function  to  check  the  bounds  fulfilment  at  the  frequencies  of 
%  operation. 

o, 

0 

%  Authors:  Mario  Garcia-Sanz,  Carlos  Molins 
%  1 . September . 2 00 9 


function  [v,  war]  =  BValid  (kg,  cz  ,  r  z ,  pc,  pr ,  bdb_a,  bdb_b,  w,  nPO ,  dPO  ) 

%  -  Bounds  restrictions  - 

Vw  =  [ ]  ; 

for  k  =  1:  length (w) 

[mod,  ph]  =  TLo (kg, cz , r z , pc, pr , nPO , dPO , w (k) ) ; 
mod  =  loglO (mod) ; 

qu  =  ilingrad ( logl 0 (bdb_a ( : , k) ) ,  ph) ; 
ql  =  ilingrad ( logl 0 (bdb_b (:, k) ) ,  ph) ; 
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Vu  =  max(qu  -  mod,  0 )  ; 

VI  =  max (mod  -  ql,  0 ) ; 
if  qu  >=  ql 

Vw  =  [Vw,  min(Vu,  VI)]; 
else 

Vw  =  [Vw,  max(Vu,  VI)]; 
end; 
end; 

%  -  Sum  and  variance  calculation 

V  =  0;  %  Bounds  violation 

war  =  0;  %  Violation  variance 

V  =  sum (Vw . ^2 ) ; 
if  length (Vw)  >  2 

war  =  var  (Vw)  ; 

end; 


%  CONFIGURATION  Configuration  of  all  the  parameters  needed  for  the 
%  automatic  loop-shaping. 

o, 

0 

%  Authors:  Mario  Garcia-Sanz,  Carlos  Molins 
%  1 . September . 2 00 9 


Starting  parameters  for  G(s) 


nrz  =  2; 
ncz  =  0; 
npr  =  1; 
npc  =  1 ; 

o 

%  Numbe 
%  Numbe 
%  Numbe 
%  Numbe 

0 

al  =  0.001; 

o, 

0 

a2  =  0.001; 

o, 

0 

a  3  =  2  ; 

o, 

0 

rl  =  [5, 

18]  ; 

o, 

0 

r2  =  [4, 

16]  ; 

o, 

0 

r3  =  [4, 

16]  ; 

o, 

0 

o 

0 

-  E- 

lambda  = 

0.2; 

o, 

0 

prob_mut 

=  0. 

08; 

prob_cros 

=  0 

.25 

-  Multi-objective  function  coefficients 

%  Gain  objective. 

%  Closeness  to  UHFB  objective. 

%  Decreasing  phase  objective. 

%  Bounds  non-fulfillment  penalty. 

%  Instability  penalty. 

%  Non-decreasing  |Lo(s) |  penalty. 


%  Mutation  probability. 
%  Cross  probability. 


%  -  Number  of  individuals  in  the  population  - 

N_vectors  =  80;  %  Number  of  individuals  in  each  generation. 

N_generations  =  300;  %  Number  of  generations. 


%  -  Poles  and  zeros  ranges 


%  Zeros  range. 
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Range_z  =  [0.01  0.01  ;  ... 

400  400  ] ; 

%  Poles  range. 

Range_p  =  [0.01  0.01  0.01  ;  ... 

1500  1500  1500  ]; 

%  Gain  range. 

Range_k  =  [0.01  1500]; 

% - Range  union 


Range  =  [Range_k',  Range_z,  Range_p] ; 


o, 

0 


o, 

0 

o, 

0 

o, 

0 


CONTROLLERS  Program  to  store  all  the  controllers  that  have  been 

Authors:  Mario  Garcia-Sanz,  Carlos  Molins 
1 . September .2009 


obtained . 


for  h=l:l:z-l  %  for  the  z-1  controllers  that  have  been  obtained. 
mrealpoles=l ;  %  Initialize  parameters 

mrealzeros=l ; 
mcomplexpoles=l ; 
mcomplexzeros=l ; 

for  hz=l : 1 :numel (find (realzeros (h, : ) >0) )  %  for  the  nzr  real  zeros 

if  realzeros (h, hz ) ~=0 

mrealzeros=conv (mrealzeros, [1/realzeros (h, hz)  1]); 

end 

end 

for  hp=l : 1 : numel ( f ind ( realpoles (h, : ) >0 ) )  %  for  the  npr  real  poles 

if  realpoles (h, hp) ~=0 

mrealpoles=conv (mrealpoles , [1/realpoles (h, hp)  1]); 

end 

end 

for  hcz=l : 2 : numel ( find (complexzeros (h, :) >0 )) -1  %for  the  complex  zeros 
if  complexzeros (h, hcz+1 ) ~=0 
omeganl=sqrt  (  (complexzeros  (h,hcz)  ''2)  + .  .  . 

(complexzeros  (h, hcz  +  1)  ''2)  )  ; 
mcomplexzeros=conv  (mcomplexzeros ,  [1/  (omeganl''2)  .  .  . 

2  *  complexzeros  (h,  hcz )  /  (omeganl  ''2 )  1  ]  )  ; 

end 

end 

for  hpc=l : 2 : numel ( find (complexpoles (h, :) >0 )) -1  %for  the  complex  poles 
if  complexpoles (h, hpc+1 ) ~=0 
omegan2=sqrt  (  (complexpoles  (h,  hpc)  ''2)  +  .  .  . 

(complexpoles  (h,  hpc+1 )  ''2  )  )  ; 
mcomplexpoles=conv  (mcomplexpoles ,  [1/  (omegan2''2)  .  .  . 

2  *  complexpoles  (h,  hpc)  /  (omegan2  ''2 )  1  ]  )  ; 

end 

end 

numG=conv (gain (h) , conv (mrealzeros, mcomplexzeros) )  ; 
denG=conv (mrealpoles , mcomplexpoles )  ; 

G (h) =tf (numG, denG) ;  %  To  store  the  z-1  controllers 


end 
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o, 

0 


o, 

0 

o, 

o 

o, 

0 

o, 

0 


COSTS  Cost  function  for  the  G(s)  controller,  considering  complex  poles 
and  zeros.  Function  to  evaluate  the  cost  function. 

Authors:  Mario  Garcia-Sanz,  Carlos  Molins 
1 . September .2009 


function  [y,  v]  =  costs (X, nr z , ncz , npr, npc, bdb_a, bdb_b, w, nPO , dPO ) 

%  -  Poles  and  zeros  extraction  from  X  vector  - 

%  Structure: 

%  [gain,  complex  zeros,  real  zeros, 

%  complex  poles,  real  poles] 


kg  =  X ( 1 )  ; 

cz  =  X (2 : l+2*ncz) ; 

rz  =  X (2+2*ncz : l+2*ncz+nrz)  ; 

pc  =  X (nr z+2 *ncz+2 : 1+nr z+2 *ncz  +  2 *npc)  ; 

pr  =  X (2+nrz+2*ncz+2*npc : l+nrz+2*ncz+2*npc+npr) ; 

%  Bounds  fulfilment 

V  =  BValid (kg, cz , r z , pc, pr , bdb_a, bdb_b, w, nPO , dPO ) ; 

%  System  stability 

St  =  stable (kg, cz , rz , pc, pr, nPO , dPO ) ; 

%  Controller  with  decreasing  |Lo(s) | 

ds  =  chkDescen (kg,  cz,rz,pc,pr,  nPO,  dPO,  [0.1  1000]); 
%  Controller  tangent  to  UHFB 

mf  =  LUHFB (kg, cz , rz , pc, pr , bdb_a, bdb_b, w, nPO , dPO , 1 00 )  ; 

%  Controller  with  decreasing  phase 

dsph  =  chkPhDes (kg,  cz,rz,pc,pr,  nPO,  dPO,  w)  ; 

%  Cost  vector 

y  =  [kg^2,  mf,  dsph,  v,  st,  ds]; 


%  CHKDESCEN  To  check  if  a  controller  produces  a  decreasing  |Lo(s) | . 

o, 

o 

%  Authors:  Mario  Garcia-Sanz,  Carlos  Molins 
%  1 . September . 2 00 9 


function  y  =  chkDescen (kg,  cz,  rz,  pc,  pr,  nPO,  dPO,  w) 

%  Open  loop  function  Lo(s) 

[numG,  denG]  =  TG (kg, cz , r z , pc, pr ) ; 

T2  =  conv(numG,  nPO); 

T1  =  conv(denG,  dPO); 

P  =  tf (T2, Tl) ; 

Mod  =  20*logl0 (abs (P) ) ; 
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wdB  =  loglO (freq) ; 
for  i=2 : length (wdB) 

k(i-l)  =  (Mod(i)  -  Mod ( i-1 ) ) / { wdB { i )  -  wdB{i-l)); 

end; 

y  =  0; 

if  max (k)  >  0 

y  =  1; 

end; 


%  CHKPHDES  Function  to  check  if  the  phase  is  always  decreasing. 

o, 

0 

%  Authors:  Mario  Garcia-Sanz,  Carlos  Molins 
%  1 . September . 2 00 9 


function  y  =  chkPhDes (kg,  cz,  rz,  pc,  pr,  nPO,  dPO,  w) 

%  Lo(s)  transfer  function 

[numG,  denG]  =  TG (kg, cz , r z , pc, pr ) ; 

T2  =  conv(numG,  nPO); 

T1  =  conv(denG,  dPO); 

%  Check  if  the  phase  increases  in  any  point. 

P  =  tf (T2, Tl) ; 

Ph  =  angle  (P ) ; 
for  i=2 : length (w) 

k (i-1)  =  Ph (i)  -  Ph (i  -  1); 
if  k (i  -  1)  <0 
k(i  -  1)  =  0; 

end; 

end; 

y  =  sum (k) ; 


%  DISORDER  Function  to  disorder  the  rows  of  a  matrix. 

o, 

0 

%  Authors:  Mario  Garcia-Sanz,  Carlos  Molins 
%  1 . September . 2 00 9 


function  Y  =  disorder (X) 


Y  =  [  ]  ; 

[r,  c]  =  size  (X)  ; 
if  r  >  0 

for  i=l : r 


[rr,  cc] 

=  size (Y) ; 

if  rr  > 

0 

pos 

=  round ( (rr 

else 

pos 

=  1; 

end; 
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Y  =  [Y (1 :pos-l, : ) ;  X{i,:);  Y (pos : rr , : ) ] ; 

end; 

end; 


%  DOMINANCE  Function  to  know  the  dominance  of  the  objectives  or 
%  constraints  in  the  fitness  function. 

o, 

0 

%  Authors:  Mario  Garcia-Sanz,  Carlos  Molins 
%  1 . September . 2009 


function  [dtl,  dt2]  =  dominance (XCost , al , a2 , a3 , rl , r2 , r3 ) 


%  Cost 

%  [kg,  Phase_dif,  Des_phase,  Bounds,  Stability,  Des_mod] 


[r,c]  =  size(XCost); 


parameters  =  [al, 
for  i=l:r 

dr 1  =  0 ; 
dr2  =  0; 
dr 3  =  0; 
if  XCost (1,4) 
dr 1  =  1 ; 

end; 

if  XCost (1,5) 
dr2  =  1; 

end; 

if  XCost (1,6) 
dr 3  =  1; 

end; 


a2,  a3,  rl,  r2, 

>  0 

>  0 

>  0 


r3]  ; 


1,  0]; 

dr3 ] . *XCost_t ( i , 4 : 6 ) ] ; 


F  =  sum (XCost_t ( : , 4 : 6 )  ,  2 )  ; 
Q  =  sum (XCost_t ( : , 1 : 3 )  ,  2 )  ; 

mu  =  F . / Q ; 

n  f  =  0  ; 
nq  =  0  ; 

for  i=l : length (mu) 


i f  mu ( i 

)  >  1 

nf 

=  n  f  +  1  ; 

else 

nq 

=  nq  +  1; 

end; 

end; 

dt2 

=  mu ( 1 ) 

r 

dtl 

=  [nq. 

nf  ]  ; 

XCost_t(i,:)  =  XCost(i,:).*parameters; 
XCost_t(i,:)  =  XCost_t(i,:)  +  [0,  0,  0,  1, 
XCost_t(i,:)  =  [XCost_t ( i , 1 : 3 ) ,  [drl,  dr2, 

end; 
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%  EvAl  Matlab  script  to  study  the  evolutive  strategies. 

o, 

o 

%  Authors:  Mario  Garcia-Sanz,  Carlos  Molins 
%  1 . September . 2 00 9 


Initial  population 


X  =  [  ]  ; 

for  i=l : N_vectors 
X_r  =  [ ]  ; 

for  j=l : length (Range) 

X_r  =  [X_r,  Range (l,j)  +  rand { 1 )* (Range { 2 , j )  -  Range ( 1 , j ))] ; 

end; 

X  =  [X;  X_r] ;  %  matrix  that  stores  in  each  row  an  individual  of  the 
%  initial  population. 

end; 

%  -  Bounds  into  non-logarithm  scale  - 

[r,  c]  =  size (bdb ( 1 : length (bdb) -2 ,:))  ; 
bdb_a  =  lO.-' (bdb(l:r/2,  :)  720)  ; 
bdb_b  =  lO.-' (bdb(r/2  +  l:r,  :)  720)  ; 

w  =  bdb ( length (bdb) -1 ,:) ;  %  bounds  frequency  vector 

%  -  Nominal  Plant  - 

nompt  =1;  %  The  nominal  plant  is  always  in  the  first  row. 

nPO  =  np ( nompt, :); 
dPO  =  dp ( nompt,:); 

%  -  Evolutive  algorithm  - 

[best_vector ,  best_ac,  dist_ac,  X,  dominant]  =  ... 

evolStr (X, N_generations , prob_mut,  prob_cros, Range, nrz, ncz, npr, npc, . . . 
al,a2,a3,rl,r2,r3, lambda, bdb_a, bdb_b, w, nPO , dPO ) ; 

%  -  Comparisons  - 

y  =  costs (best_vector, nrz , ncz , npr , npc, bdb_a, bdb_b, w, nPO , dPO ) ; 

message (' Best  controller  found  (X) Cost :', best_vector, nrz , ncz , npr , npc, y) ; 

%  -  Controller  validity  - 


kg  =  best_vector ( 1 )  ; 

cz  =  best_vector (2 : l+2*ncz)  ; 

rz  =  best_vector (2+2*ncz : l  +  2*ncz+nrz)  ; 

pc  =  best_vector (nrz+2*ncz+2 : l+nrz+2*ncz+2*npc) ; 

pr  =  best_vector (2+nrz+2*ncz  +  2*npc : l+nrz  +  2*ncz  +  2*npc+npr)  ; 

St  =  stable (kg, cz , rz , pc, pr, nPO , dPO )  ; 

if  St  >  0 

disp (' unstable  controller'); 

else 

disp ('stable  controller'); 

end 
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%  EVOLSTR  Function  to  implement  a  simple  evolutive  strategy. 

o, 

0 

%  Authors:  Mario  Garcia-Sanz,  Carlos  Molins 
%  1 . September . 2009 


function  [best_vector ,  best_ac,  dist_ac,  X,  dominant]  =  ... 

evolStr (X, M_count, prob_mut,  prob_cros.  Range, nrz, ncz, npr, npc, . . . 
al,a2,a3,arl,ar2,ar3, lambda, bdb_a, bdb_b, w, nPO ,  dPO ) 

rl  =  arl  ( 1 )  ; 
r2  =  ar2  ( 1 ) ; 
r3  =  ar3  ( 1 )  ; 

rand ( ' state ' , 0 ) ; 

N_vectors  =  size (X ( : , 1 ) ) ; 

%  -  Evolutionary  Estrategy  (ES)  - 

disp(' Start  of  the  evolutive  method  (ES)  ...  ' ) ; 

disp  (  '  '  ) ; 

%  -  Acumulators  setup  - 

dist_ac  =  [  ] ; 
best_vector  =  [  ] ; 
best_ac  =  [ ] ; 
mut_number  =  0 ; 
mut_success  =  0; 
dominant  =  [  ] ; 

%  -  Best  vector  of  the  initial  population  - 

F  =  [  ]  ; 

D  =  [  ]  ; 

for  i=l : N_vectors 

[F(i,:),  D(i)]  =  costs (X ( i ,:), nrz , ncz , npr , npc, bdb_a, bdb_b, w, nPO , dPO )  ; 
end; 

[X,  indexF,  X_cost2 ]=selection{X,F,al,a2,a3,rl,r2,r3,rl,r2,r3,0, N_vectors ) ; 

best_vector  =  X(l,:); 

best_ac  =  [best_ac;  best_vector] ; 

dist_ac  =  X_cost2 ( 1 , 4 ) ; 

score  =  indexF ( 1 , : ) ; 

[dtl,  dt2]  =  dominance (X_cost2 , al , a2 , a3 , rl ,  r2  ,  r3 )  ; 
dominant  =  [dominant;  [dtl,  dt2]]; 

%  -  Operation  cycle  - 

rp  =  [arl;  ar2;  ar3]; 

count  =  0;  %  Iterations  accountant 

while  (count  <  M_count) 

%  -  Penalties  establishment  - 
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r_new  =  rp ( : , 1 )  +  (rp(:,2)  -  rp { : , 1 ) ) *count/M_count ; 
nr 1  =  r_new ( 1 ) ; 
nr 2  =  r_new (2 ) ; 
nr 3  =  r_new ( 3 ) ; 

-  Crossover  - 

X_cost  =  [ ] ; 

X_cross  =  [];  %  Population  for  the  crossover. 

X_no_cross  =  []; 
for  i=l : N_vectors 

if  rand(l)  <  prob_cros 

X_cross  =  [X_cross;  X{i,:)]; 
else 

X_no_cross  =  [X_no_cross;  X{i,:)]; 
end; 
end; 

mark  =0;  %  If  mark  ==  1,  we  add  one  element  to  have  an  even  number. 

[r,c]  =  size(X_cross); 

1  f  r  >  0 

if  rem ( length (X_cross (:, 1 )) ,  2)  >  0 

X_cross  =  [X_cross;  X_cross { length {X_cross {:,  1 )),:)]  ; 
mark  =  1; 
end; 

X_cross=disorder (X_cross ) ; %  Disorder  the  rows  of  the  matrix. 
X_cross_new  =  []; 
for  1=1 : 2 : length (X_cross (:,  1 )) -1 
cross_element_l  =  [];  %  Child  1 
cross_element_2  =  [];  %  Child  2 
for  k=l : length (X ( 1 ,:) ) 
if  rand ( 1 )  <  0.3 
a  =  rand ( 1 ) ; 
else 

a  =  round ( rand ( 1 )) ; 
end 

cross_element_l  =  [ cross_element_l ,  a*X_cross (1, k)  +  ... 

(1-a) *X_cross (1  +  1,  k)  ]  ; 

cross_element_2  =  [ cross_element_2 ,  (l-a) *X_cross (1, k)  +  ... 

a  *  X_c  r o  s  s ( 1  + 1 ,  k )  ]  ; 

end; 

%  Join  the  parents  with  their  children  and  select  the  two  bests. 
X_sub_cross  =  [X_cross { 1 , : ) ;  X_cross { 1+1 , : ) ;  cross_element_l;  ... 
cross_element_2 ] ; 

F  =  [  ]  ; 

D  =  [  ]  ; 
for  1=1:4 

[F(i,:)  D(i)]  =  costs {X_sub_cross { 1 ,:), nrz , ncz , npr, npc, bdb_a, .. . 
bdb_b , w , nP  0 , dP  0 )  ; 

end; 

[X_cross2,  I_cross,  X_cost2]  =  selection {X_sub_cross,  F,  al,a2,... 

a3, rl, r2, r3,nrl,nr2,nr3, 2, 2)  ; 

X_cross_new  =  [X_cross_new;  X_cross2]; 

X_cost  =  [X_cost;  X_cost2]; 
end; 
end; 

%  -  Mutation  - 
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X_mut  =  [  ] ; 

X_no_mut  =  [ ] ; 

[r,  c]  =  size  (X)  ; 
for  i  =  1 : r 

if  rand(l)  <  c*prob_mut 
if  rand ( 1 )  >  0.3 

%  Only  one  parameter  mutates . 

pos  =  round ( (c-1 ) *rand { 1 ) +1 ) ;  %  Mutation  position. 

new_value  =  randomr (X { i , pos ) ,  Range { 1 ,  pos ),  Range ( 2  ,  pos )  ,  lambda); 

Mut_vector  =  X(i,:); 

Mut_vector (pos )  =  new_value; 

else 

%  All  the  parameters  mutate, 
for  k  =  1 : c 

Mut_vector (k)  =  randomr (X (i, k) , Range (1, k) , Range (2, k) , lambda) ; 

end; 

end; 

F  =  [  ]  ; 

D  =  [  ]  ; 

[F(l,  :) ,D(1) ]=costS (Mut_vector , nrz,ncz,npr,npc, bdb_a, bdb_b, w, nPO , dPO )  ; 
[F(2,  :) ,D(2) ]=costS(X{i,  :) ,nrz,ncz,npr,npc, bdb_a, bdb_b, w, nPO , dPO ) ; 
[X_mut2 , I_mut , X_cost2 ] =selection { [Mut_vector ;  X{i,  :)  ],F,al,a2,a3,  .  .  . 

rl,r2,r3,nrl,nr2,nr3,l,l) ; 

X_mut  =  [X_mut;  X_mut2]; 

X_cost  =  [X_cost;  X_cost2]; 

if  X_mut2  ==  Mut_vector 

mut_success  =  mut_success  +  1; 

end; 

else 

X_no_mut  =  [X_no_mut;  X{i, :) ] ; 
end; 
end; 

[r,  c]  =  size(X_mut); 
mut_number  =  mut_number  +  r; 

%  -  Evolution  of  mutation  range  - 

if  rem(count,  5)  ==  0 

disp ( sprintf ( ' Step  number  %d.  Mutation  success  %2 . 3f count, .. . 

mut_success/mut_number) )  ; 
vari_s  =  [ ' ' ] ; 
param_s  =  [ ' ' ] ; 
score_s  =  [ ' ' ] ; 
dominant_s  =  [  '  ' ] ; 
for  i=l : length (best_vector) 

param_s  =  sprintf ('%s  %4.4f',  param_s,  best_vector ( i ) ) ; 
vari_s  =  sprintf ('%s  %2.4f',  vari_s,  sqrt (var (X ( : , i ) ) ) ) ; 

[sr,  sc]  =  size (score) ; 

end; 

for  k=l : sc 

score_s  =  sprintf ('%s  %4.3f  score_s,  score (sr, k) ) ; 

end; 

[dr,dc]  =  size (dominant) ; 
for  k=l : dc 

dominant_s  =  sprintf ('nq  =  %d,  nf  =  %d,  %3 . 4f ', dominant (dr , 1 ),.. . 

dominant (dr, 2) ,  dominant (dr, 3) ) ; 
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end; 

disp ( sprintf (' Vector :  %s',  param_s)); 

disp ( sprintf (' Stand . Dev . :  %s',  vari_s)); 
disp (sprintf ('Classification:  [%s] score_s ) ) ; 
disp ( sprintf (' Population  balance.:  %s',  dominant_s ) ) ; 
disp  (  '  '  ) ; 

mut_success  =  0; 
mut_number  =  0; 
end; 

%  -  Selection  - 

%  X  =  [X_no_mut;  X_cross_new;  X_mut ] . 

X  =  [X_cross_new;  X_mut]; 
nn  =  size (X ( : , 1 ) ) ; 

X  =  [X;  X_no_mut]; 

[mr,  me]  =  size (X_no_mut) ; 

F  =  [  ]  ; 
for  i=l :mr 

F(i,:)  =  costs (X_no_mut (i,  :) ,nrz,ncz,npr,npc, bdb_a, bdb_b, w, nPO , dPO )  ; 
end; 

%  F  =  [F;  X_cost] . 

F  =  [X_cost;  F]; 

[X_new,  I_new,  Cost_new]  =  selection (X, F, al, a2, a3, rl, r2, r3, nrl, nr2, .. . 
nr3,  nn, N_vectors) ; 

%  Population  update. 

X  =  X_new; 

%  Acumulators  update. 

best_vector  =  X(l,:); 

dist_ac  =  [dist_ac,  Cost_new ( 1 , 4 ) ] ; 

score  =  [score;  I_new(l,:)]; 

best_ac  =  [best_ac;  best_vector] ; 

[dtl,  dt2]  =  dominance (Cost_new,  al , a2 , a3 , nrl , nr2 , nr3 ) ; 
dominant  =  [dominant;  [dtl,  dt2]]; 

%  Penalties  update, 
r  1  =  n  r  1  ; 
r2  =  nr2; 
r  3  =  n  r  3  ; 

%  Accountant  increment, 
count  =  count  +  1; 

end; 

%  ILINGRAD  Function  for  the  linear  interpolation  of  a  function  (Fin)  with 
%  repsect  to  the  phase  from  0  to  -360  degrees,  in  a  x  point. 

o, 

0 

%  Authors:  Mario  Garcia-Sanz,  Carlos  Molins 
%  1 . September . 2009 

o, 

0 
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function  result  =  ilingrad (Fin,  x) 

%  Number  of  points  of  Fin 
N  =  length (Fin)  ; 

%  Function  vector  invertion 
for  i=l : N 

F (i)  =  Fin (N+l-i) ; 
end; 

%  Phase  division  in  degrees, 
delta  =  -360/ (N-1); 
na  =  N  -  f loor (x/delta) ; 
nb  =  N  -  ceil (x/delta) ; 

%  Linear  interpolation 
if  na  ==  nb 

result  =  F (nb) ; 
else 

result  =  F (nb)  +  (N  -  x/delta  -  nb) * (F (na)  -  F (nb) ) ; 
end; 

%  Maximum  and  minimum  value  limitation, 
if  result  >  0 

if  result  <  0.00001 
result  =  0.00001; 
end; 
else 

if  result  >  -0.00001 
result  =  -0.00001; 
end; 
end; 


%  ILINMOD  Function  for  the  linear  interpolation  of  a  function  (Fin)  with 
%  respect  to  the  module. 

o, 

0 

%  Authors:  Mario  Garcia-Sanz,  Carlos  Molins 
%  1 . September . 2 00 9 


function  result  =  ilinmod(Fin,  x) 

%  Number  of  points  of  Fin 
N  =  length (Fin  ( 1 ,:))  ; 

Fmod  =  Fin  ( 1 ,  :  ) ; 

FPh  =  -  Fin  (2,  :  ) ; 

if  (x  <  min (Fmod) )  |  (x  >  max (Fmod)) 

result  =  -1; 

else 

i  =  1; 

while  X  >=  Fmod(i) 
i  =  i  +  1  ; 

end; 

ml  =  Fmod ( i  -  1 ) ; 
phi  =  FPh (i  -  1)  ; 
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i  =  N; 

while  X  <=  Fmod(i) 
i  =  i  -  1; 

end; 

m2  =  Fmod ( i  +  1 )  ; 
ph2  =  FPh ( i  +  1); 

if  ml  <  m2 

result  =  phi  +  (x  -  ml)  *  {ph2  -  phi)  /  {m2  -ml)  ; 

else 

result  =  phi; 

end; 

end; 


%  K_COR  Function  to  correct  the  value  of  Kg  due  to  the  way  the  poles  and 
%  zeros  are  defined. 

o, 

0 

%  Authors:  Mario  Garcia-Sanz,  Carlos  Molins 
%  1 . September . 2 00 9 


function  Kg  =  K_cor (kg, cz , r z , pc, pr ) 

Kg  =  kg; 

numG  =  kg; 
denG  =  1.0; 

%  Complex  zeros 

for  j=l : 2 : numel ( f ind (cz ( : ) >0 ) ) 

Kg  =  Kg  /  (cz(j)^2  +  cz{j+l)^2); 

end 


%  Real  zeros 

for  j  =  l : numel ( find ( rz (:) >0  )  ) 

Kg  =  Kg  /  rz  (  j ) ; 

end; 

%  Complex  poles 

for  j  =  l : 2 : numel ( f ind (pc  (  : ) >0 ) ) 

Kg  =  Kg  *  (pc(j)^2  +  pc{j+l)^2); 

end; 

%  Real  poles 

for  j=l : numel ( find (pr (:) >0 ) ) 

Kg  =  Kg  *  pr  (  j ) ; 

end; 
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%  K_COR2  Function  to  correct  the  value  of  Kg  due  to  the  way  the  poles  and 
%  zeros  are  defined. 

o, 

o 

%  Authors:  Mario  Garcia-Sanz,  Carlos  Molins 
%  1 . September . 2 00 9 


function  Kg  =  K_cor2 (kg, cz , r z , pc, pr ) 

Kg  =  kg; 

numG  =  kg; 
denG  =  1.0; 

%  Complex  zeros 

for  j=l : 2 : numel ( f ind (cz ( : ) >0 ) ) 

Kg  =  Kg  *  (cz(j)^2  +  cz(j+l)^2); 

end 


%  Real  zeros 

for  j=l :numel (find (rz ( : ) >0) ) 

Kg  =  Kg  *  r  z  (  j ) ; 

end; 

%  Complex  poles 

for  j  =  l : 2 : numel ( f ind (pc  (  : ) >0 ) ) 

Kg  =  Kg  /  (pc(j)^2  +  pc(j  +  l)''2); 

end; 


%  Real  poles 

for  j=l : numel ( find (pr (:) >0 ) ) 

Kg  =  Kg  /  pr  (  j ) ; 

end; 

%  LUHFB  Function  to  evaluate  how  tangent  is  Lo  to  the  UHFB . 

o, 

0 

%  Authors:  Mario  Garcia-Sanz,  Carlos  Molins 
%  1 . September . 2009 

o, 

0 


function  ph_ac  =  LUHFB (kg, cz , rz , pc, pr , bdb_a, bdb_b,  w,  nPO ,  dPO ,  w_ref ) 

%  -  Bounds  restrictions  - 

ph_ac  =  0; 

for  k  =  1: length (w) 

[mod,  ph]  =  TLo (kg, cz , r z , pc, pr , nPO , dPO , w (k) ) ; 
mark  =  min (20*logl0 (bdb_b ( : , k) ) ) ; 
mod  =  20*logl0 (mod) ; 
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if  (mark  <  320)  &  (w(k)  >=  w_ref) 

Y  =  uph (20*logl0 (bdb_a) ,  20*logl0 {bdb_b) ) ; 
if  mod  >  min(y(l,:)) 

ph_b  =  ilinmod(y,  mod)  ; 
if  ph_b  >  0 

diff  =  abs (ph_b  -  abs (ph) ) ; 
if  diff  >  5 

ph_ac  =  ph_ac  +  diff; 

end; 

end; 

end; 

end; 

end; 


%  MENUPLOT  Menu  to  choose  the  task  to  do  with  the  results  obtained  from  the 
%  automatic  loop-shaping. 

o, 

0 

%  Authors:  Mario  Garcia-Sanz,  Carlos  Molins 
%  1 . September . 2009 


k=menu ( ' Choose  an  option ' ,  ' step ' ,  ' bode ' ,  ' nichols ' ,  ' nyquist ' ,  ' chksiso ' )  ; 

switch  k 

case  1 

for  h=l : 1 : z-1 
figure  (2 ) 

step  (cp  (1,  1,  1)  *G  (h)  )  ; 
grid  on; 

set (gca,  ' FontSize ' ,  [5] ) 
axis  (  [0  70  0  200] ) 
hold  all 

end 
case  2 

for  h=l : 1 : z-1 
figure  (2 ) 

bode (cp ( 1 , 1 , 1 ) *G (h) , w) ; 
grid  on; 

set (gca,  ' FontSize [5] ) 
axis([0.001  10000  -180  -45]) 
hold  all 

end 
case  3 

load  wloop; 
for  h=l : 1 : z-1 

Ipshape (wloop,  bdb,  cp (1,1,1),  G(h),  []); 

end 
case  4 

for  h=l : 1 : z-1 
figure  (2 ) 

nyquist (cp ( 1 , 1 , 1 ) *G (h) ) ; 
grid  on; 

set (gca,  ' FontSize ' ,  [5] ) 
axis ( [-6  6  -3  3 ] ) 
hold  all 

end 
case  5 

load  wcheck; 
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load  W1 
load  W7 
for  h=l : 1 : z-1 

If  ~isempty (W1 ) 

chksiso ( 1 , wcheck, Wl,cp,0,G{h),l,l); 

title ([' bound  type  1  with  the  ' , num2str (h) , ' th  controller']); 

grid  on; 

end 

if  ~isempty (W2 ) 

chksiso (2 , wcheck, W2,cp,0,G{h),l,l); 

title ([' bound  type  2  with  the  ' , num2str (h) , ' th  controller']); 

grid  on; 

end 

if  -isempty (W3 ) 

chksiso ( 3 , wcheck, W3,cp,0,G{h),l,l); 

title ([' bound  type  3  with  the  ' , num2str (h) , ' th  controller']); 

grid  on; 

end 

if  ~isempty (W4 ) 

chksiso ( 4 , wcheck, W4,cp,0,G{h),l,l); 

title ([' bound  type  4  with  the  ' , num2str (h) , ' th  controller']); 

grid  on; 

end 

if  ~isempty (W5 ) 

chksiso ( 5 , wcheck, W5,cp,0,G{h),l,l); 

title ([' bound  type  5  with  the  ' , num2str (h) , ' th  controller']); 

grid  on; 

end 

if  ~isempty (W6 ) 

chksiso ( 6 , wcheck, W6,cp,0,G{h),l,l); 

title ([' bound  type  6  with  the  ' , num2str (h) , ' th  controller']); 

grid  on; 

end 

if  ~isempty (W7 ) 

chksiso ( 7 , wcheck, W7,cp,0,G{h),l,l); 

title ([' bound  type  7  with  the  ' , num2str (h) , ' th  controller']); 

grid  on; 

end 

if  ~isempty (W8 ) 

chksiso ( 8 , wcheck, W8,cp,0,G{h),l,l); 

title ([' bound  type  8  with  the  ' , num2str (h) , ' th  controller']); 

grid  on; 

end 

if  ~isempty (W9 ) 

chksiso ( 9 , wcheck, W9,cp,0,G{h),l,l); 

title ([' bound  type  9  with  the  ' , num2str (h) , ' th  controller']); 

grid  on; 

end 

end 

end 


%  MESSAGE  Function  to  show  the  gain, 

o, 

0 

%  Authors:  Mario  Garcia-Sanz,  Carlos 
%  1 . September . 2009 


the  zeros,  the  poles  and  the  cost 
Molins 


function  result  =  message (ml, m2, x, nrz, ncz, npr, npc, cost) 
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disp (ml ) ; 

kg  =  X ( 1 ) ; 

cz  =  X (2 : l+2*ncz) ; 

rz  =  X (2+2*ncz : l+2*ncz+nrz)  ; 

pc  =  X (nrz+2*ncz+2 : l+nrz+2*ncz  +  2*npc)  ; 

pr  =  X (2+nrz+2*ncz+2*npc : l+nrz+2*ncz+2*npc+npr) ; 

khf  =  K_cor(kg,  cz,  rz,  pc,  pr)  ; 

S  =  sprintf ( '  kg  =  %4.2f  <HFG  =  %4.2f>\n  Zeros:  ',  kg,  khf )  ; 
for  i=l : size (rz, 2) 

S  =  [S,  sprintf (' %3 . 2f  ',rz(i))]; 
end; 

for  i=l : 2 : size (cz, 2) 

omega  =  sqrt(cz(i)^2  +  cz  ( i  +  1 ) ''2  )  ; 
chi  =  cz (i) /omega; 

S  =  [S,  sprintf (' (%3.2f  %3.2f)  [wn  =  %3.2f,  chi  =  %3.2f]  ',cz(i), 

cz(i+l),  omega,  chi) ] ; 
end; 

S  =  [S,  sprintf (' \n  Poles;  ')]; 
for  i=l : size (pr, 2) 

S  =  [S,  sprintf (' %3 . 2f  ',pr(i))]; 
end; 

for  i=l : 2 : size (pc, 2) 

omega  =  sqrt(pc(i)^2  +  pc (i+1) ^2); 
chi  =  pc (i) /omega; 

S  =  [S,  sprintf (' (%3.2f  %3.2f)  [wn  =  %3.2f,  chi  =  %3.2f]  ',pc(i), 

pc (i+1),  omega,  chi)]; 
end; 

disp ( S ) ; 
result  =  S; 

S  =  sprintf  (  '  %s',  m2 )  ; 

for  i=l : length (cost) 

S  =  [S,  sprintf (' %4 . 5f  ',  cost(i))]; 

end; 

disp ( S ) ; 

result  =  [result,  S]  ; 

%  ORDER4  Indexes  weighting  based  on  different  criteria. 

o, 

0 

%  Authors:  Mario  Garcia-Sanz,  Carlos  Molins 
%  1 . September . 2 00 9 

o, 

0 


function  [y,  IM]  =  order4  (X,  al,  a2,  a3,  rl,  r2,  r3,  nrl,  nr2,  nr3,  nn)  ; 

[r,  c]  =  size  (X)  ; 

%  Indexes  weighting 
for  k=l : r 

drl  =  0; 
dr2  =  0; 
dr 3  =  0; 

if  X(k,4)  >0  %  Bounds  violation 

drl  =  1; 
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end; 

if  X(k,5)  >0  %  Stability  violation 

dr2  =  1 ; 

end; 

if  X(k,6)  >0  %  Lo  with  decreasing  module  violation 

dr3  =  1; 

end; 

if  k  >  nn 

IM(k,:)  =  [X (k, 1 : 3) . * [al  a2  a3],  [drl  dr2  dr3].*[rl  r2  r3].*... 
([1/rl  l/r2  l/r3] .*X{k,4:c)  +  1)]; 

else 

IM(k,:)  =  [X (k, 1 : 3) . * [al  a2  a3],  [drl  dr2  dr3].*[nrl  nr2  nr3].*... 
([1/nrl  l/nr2  l/nr3] . *X (k, 4 : c)  +  1)]; 

end; 

end; 


MIV  =  sum(IM,  2);  %  Sum  the  elements  of  each  row. 

%  To  put  in  order  and  obtain  the  final  order  index. 
[Ord2,  III]  =  sort (MIV); 


y  =  III; 

for  k=l : r 

y (III  (k) )  =  k; 

end; 


%  RANDOMR  Function  that  selects  a  random  number,  uniformly  distributed, 
%  within  the  interval . 

o, 

0 

%  Authors:  Mario  Garcia-Sanz,  Carlos  Molins 
%  1 . September . 2 00 9 


function  res  =  randomr (x,  rangl,  rang2,  lambda) 

%  Total  range 

rang  =  (rang2  -  rangl); 

%  Result  limitation  within  the  range 
res  =  X  +  lambda* ( rand ( 1 )  -  l)*rang; 
if  rangl  >  res 
res  =  rangl; 

end; 

if  rang2  <  res 
res  =  rang2; 

end; 


%  SELECTION  Function  to  select  the  individuals  of  a  population  in  terms  of 
%  their  evaluation  in  the  fitness  function. 

o, 

0 

%  Authors:  Mario  Garcia-Sanz,  Carlos  Molins 
%  1 . September . 2009 


function  [y,  ly,  Cy]  =selection  (X,  XCost,  al,  a2,  a3,  rl,  r2,  r3,  nrl,  nr2,  nr3,  nn,  n) 
[FF,  I]  =  order4  (XCost ,  al,  a2,  a3,  rl,  r2,  r3,  nrl,  nr2,  nr3,nn); 
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ly  =  []; 

y  =  [  ]  ; 

Cy  =  [ ]  ; 
for  i=l : n 

pos  =  find(FF==i); 
y  =  [y;  X (pos, : ) ] ; 
ly  =  [ly;  I (pos, : ) ] ; 

Cy  =  [Cy;  XCost (pos , : ) ] ; 

end; 


%  STABLE  Stability  of  the  1  +  G(s)Po(s)  characteristic  equation.  If  any 
%  root  lies  on  the  RHP  the  system  is  instable. 

o, 

0 

%  Authors:  Mario  Garcia-Sanz,  Carlos  Molins 
%  1 . September . 2009 


function  st  =  stable (kg, cz , rz , pc, pr, nPO , dPO ) 

%  Transfer  function  for  Lo(s) 

[numG,  denG]  =  TG (kg, cz , r z , pc, pr ) ; 

T2  =  conv(numG,  nPO); 

T1  =  conv(denG,  dPO); 

%  Characteristic  equation  calculation 
rl  =  length (T1 ) ; 
r2  =  length  (T2); 
if  (rl>r2) 

T2  =  [ zeros ( 1 , rl-r2 ) ,  T2]; 

else 

T1  =  [ zeros ( 1 , r2-rl ) ,  Tl]; 

end 

T  =  Tl  +  T2; 

%  Characteristic  equation  roots 

r=roots (T) ; 

rreal  =  max ( real ( r) ) ; 

if  rreal  >  0 

St  =  rreal; 

else 

St  =  0; 

end 

%  TG  Transfer  function  for  G(s) .The  outcome  is  the  numerator  and  the 
%  denominator. 

o, 

0 

%  Authors:  Mario  Garcia-Sanz,  Carlos  Molins 
%  1 . September . 2 0 0 9 

o, 

0 


function  [numG,  denG]  =  TG (kg, cz , r z , pc, pr ) 

%  -  Transfer  function  for  G(s)  - 

%  E  (s/z  +  1) 

%  G ( s )  =  kg  - 

%  E  (s/p  +  1) 
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numG  =  kg; 
denG  =  1.0; 


%  Complex  zeros 
for  j  =  l : 2 : size  (cz, 2) 

omega  =  cz(j)^2  +  cz(j+l)^2; 
chi  =  cz ( j )  ; 

numG  =  conv(numG,  [1 /omega,  2* chi /omega,  1]); 

end 

%  Real  zeros 

for  j  =  l : size  (rz, 2) 

numG  =  conv(numG,  [l/rz{j)  1]); 

end 


%  Complex  poles 
for  j=l : 2 : size (pc, 2 ) 

omega  =  pc(j)''2  +  pc(j  +  l)^2; 
chi  =  pc ( j ) ; 

denG  =  conv(denG,  [1/omega,  2*chi/omega,  1]); 

end 


%  Real  poles 

for  j=l : size (pr , 2 ) 

denG  =  conv(denG,  [l/pr{j)  1]); 

end 


o, 

o 


o, 

0 

o, 

o 

o, 

0 

o, 

0 

o, 

0 

o, 

o 

o, 

0 


TLo  Transfer  function  for  Lo(s)  =  G{s)Po{s) 

E  (s  +  z) 

Lo  ( s )  =  kg - Po  { s ) 

E  (s  +  p) 

Authors:  Mario  Garcia-Sanz,  Carlos  Molins 
1 . September .2009 


function  [mod,  ph]  =  TLo (kg, cz , rz , pc, pr , nPO , dPO , f req) 

[numG,  denG]  =  TG (kg, cz , rz , pc, pr) ; 
numLO  =  conv(numG,  nPO); 
denLO  =  conv(denG,  dPO); 

Imo  =  polyval (numLO ,  sqrt (-1 ) *f req)  /  polyval (denLO,  sqrt (-1 ) *f req)  ; 
mod  =  abs (Imo) ; 
ph  =  angle (Imo); 

% - Phase  between  0  and  -2pi - 

if  ph  >  0 

ph  =  -2*pi  +  ph; 
end; 

if  ph  <  -2*pi 

ph  =  ph  +  2*pi*f loor (-ph/ (2*pi) )  ; 
end; 

ph  =  ph  *  180  /  pi;  %  Convert  the  phase  into  degrees. 

%  TMPL-BNDS  Function  to  calculate  the  problem  templates  and  the  bounds. 

o, 

0 
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%  Authors:  Mario  Garcia-Sanz,  Carlos  Molins 
%  1 . September . 2 0 0 9 


function  [bdb,  np,  dp]  =  Tmpl_Bnds ( ) ; 

close  all; 

clc; 

clear; 

%  -  Templates  generation  - 

disp( 'Plant  templates  calculation  ...') 
c=l; 

for  k=[l  2.5  6  10] 

for  a=[l  1.5  3  6  10] 
np (c, : ) = [k*a] ; 
dp ( c , : ) = [ 1  a  0]; 
cp  (1,  1,  c)  =tf  (np  (c,  :  )  ,  dp  (c,  :  )  )  ; 
c=c+l ; 

end 

end 

disp('  ...  Finished'); 

%  -  Nominal  plant  parameters  - 

nompt  =1;  %  Nominal  plant  is  always  in  the  first  row. 

%  -  Frequency  range  - 

w  =[0.1  0.5  1  2  5  10  15  40  60  80  120  170  200  300  400  500  600  800  1000  ... 

1250  1500  1750  2000  3000  5000  10000]; 
plottmpl (w, cp, nompt)  ; 

set (gca,  ' Title ', title (' P (s)  plant  templates')); 
pause (2 ) ; 

% - Bounds  initialization - 


Wl=[] 
W2=[] 
W3=[] 
W4=[] 
W5=[] 
W6=  [  ] 
W7=[] 
W8=[] 
W9=[] 


% - Bound  Type  1 - 

disp('Type  1  bounds  calculation  ...'); 

Wl=l . 2 ;  save  W1 
wbdl=w; 

bdbl  =  sisobnds ( 1 , wbdl , W1 ,  cp) ; 
plotbnds (bdbl) ; 

set (gca,  ' Title ' , title ( ' Stability  bounds ' ) ) ; 

disp('  ...  Finished'); 

pause  (2 ) ; 


70 


Advanced  Quantitative  Robust  Control  Engineering:  Ref.  BOARD  Award  #  FA8655-08-1-3027 

New  Solutions  for  Automatic  Loop-shaping  for  SISO  Ref.  Grant  #  083027 

and  MIMO  Systems.Part  1:  SISO  Systems.  Date:  l.September.2009  Page  71 


% - Bound  Type  7  - 

wbdb7  =  [0.1  0.5  1  10]; 

disp('Type  7  bounds  calculation 

emu  =  tf(0.6584*[l  30],  [1  4  19.752]); 

cml  =  tf(120,  [1  17  82  120]); 

W7  =  [emu;  cml];  save  W7 

bdb7  =  sisobnds(7,  wbdb7,  W7,  cp) ; 

plotbnds (bdb7) ; 

set (gca,  ' Title ' , title ( ' Tracking  bounds ' ) ) ; 

disp('  ...  Finished'); 

pause (2 ) ; 

%  -  Bounds  intersection  - 

bdb  =  grpbnds (bdbl , bdb7 ) ; 
bdb  =  sectbnds (bdb) ; 
plotbnds (bdb) ; 

set  (gca,  ' Title title ( 'All  bounds  intersection')); 
disp( 'Bounds  calculation  finished.'); 
disp  (  '  '  ) ; 

pause  (2 ) ; 


UPH  Function  necessary  to 

Authors:  Mario  Garcia-Sanz 
1 . September .2009 


determine  how  tangent 
Carlos  Molins 


Lo  is  to  the  UHFB 


function  y  =  uph(bdb_a,  bdb_b) 

n  =  length (bdb_a ( 1 , : ) ) ; 
m  =  length (bdb_a (:,!)); 
ph  =  0:-5:-180; 

for  i=l :m 

bdba(i)  =  bdb_a (m  -  i  +  1,  n)  ; 
bdbb(i)  =  bdb_b (m  -  i  +  1,  n)  ; 

end; 


[bdb_aS,  la]  =  sort (bdba ( 1 : length (ph) )) ; 
[bdb_bS,  Ib]  =  sort (bdbb ( 1 : length (ph) )) ; 
bdb_aMin  =  bdb_aS(l); 
bdb_bMax  =  bdb_bS(l); 

for  i=l : length (ph) 

Ph_a (i)  =  ph (la (i) ) ; 

Ph_b (i)  =  ph (Ib  (i) ) ; 

end; 

bdb_aS2  =  []; 
bdb_bS2  =  []; 

Ph_a2  =  [ ]  ; 

Ph_b2  =  [  ]  ; 

for  i=l : length (bdb_aS ) 

If  abs (bdb_aS ( i ) )  ~=  320 
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bdb_aS2  =  [bdb_aS2,  bdb_aS{i)]; 
Ph_a2  =  [Ph_a2,  Ph_a{i)]; 

end; 

if  abs (bdb_bS ( i ) )  ~=  320 

bdb_bS2  =  [bdb_bS2,  bdb_bS{i)]; 
Ph_b2  =  [Ph_b2,  Ph_b{i)]; 

end; 

end; 


bdb_u  =  [bdb_bS2,  bdb_aS2]; 
Ph_u  =  [Ph_b2,  Ph_a2]; 

y  =  [bdb_u;  Ph_u] ; 
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8.  APPENDIX  B:  MATLAB,  C++  AND  TEXT  FILES  FOR  THE 
STRATEGY  2  (GENETIC  ALGORITHMS) 

8.1.  MATLAB  FILES 

%  AUTOLOOP SHAPING_AG  Program  to  sinthesys  optimal  controllers  with  an 
%  initial  structure  which  can  be  modified  by  means  of 

%  the  optimal  Hankel  norm  approximation.  If  any 

%  modification  occurs  the  program  sinthesys  the  optimal 

%  controller  for  thid  new  structure. 

o, 

0 

%  Authors:  Mario  Garcia-Sanz,  Carlos  Molins 
%  1 . September . 2009 

o, 

0 


%  Genetic  Algorithms 

%  First  of  all  the  user  has  to  configure  the  CFG. TXT,  CFGCONTR.TXT  and 
%  CFGCOST.TXT  files.  Take  into  account  that  this  program  has  been 
%  configured  for  controllers  with  up  to  five  poles  and  zeros. 


%  In  addition  the  plant  and  the  performance  specifications  have  to  be 
%  defined,  obtaining  the  plant  templates  and  the  bounds. 

Tmpl_Bnds_AG;  %  Execute  only  if  no  template  approximation  is  desired. 
LOOPSHN;  %  Execute  only  for  template  approximations. 

auxx= ' T ' ; 
auxx2=l ; 

while  auxx== ' T ' 


%  Now  the  OPTQFT.exe  is  executed.  As  a  result  the  following  files  are 
%  obtained:  AAGG.LOG,  PRIV1.TXT  and  VALIDVECTOR.TXT. 

10PTQFT.exe 


[C,B,A,D]  =  textread( 'CFGCONTR.txt' ,.. . 

'%f  %  [^'l,  2, 3, 4, 5]  %f  %  [■'1, 2, 3, 4, 5]  '  ,  '  delimiter  '  , 


B=2*A+1; 


aux2= ' T ' ; 

%  Once  the  user  has  the  PRIV1.TXT  the  results  are  stored  in  the  vectors 
%  k (gains),  zi (zeros),  pi (poles),  best_value (best  fitness  value)  and 
%  mean_value (generation  fitness  mean  value) . 

switch  B 

case  11 

[best_value, mean_value, Cl,k,zl,z2,z3,z4,z5,pl,p2,p3,p4,p5,c2]=. . . 
textread ( 'privl .txt ' , ' %f  %f  %c  %f  %f  %f  %f  %f  %f  %f  %f  %f  . . . 
%f  %n  %s delimiter ' ,  '  ' ) ; 
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[valids ] =t ext read { ' validvector . txt ' , ' %u ' ) ; 

%  valids  is  a  column  vector  with  1  for  feasible  solutions  and 
%  0  for  the  non-feasible  ones,  between  the  generations  of  this 
%  iteration. 

%  From  the  above  mentioned  vectors  only  the  feasible  solution 
%  are  of  interest,  so  the  following  lines  extract  them  and 
%  store  in  the  vectors  valid,  gain,  p  and  z. 

j=i; 

p=zeros (numel ( find (valids {:) ==1 ) ) , A) ; 
z  =  zeros (numel (find (valids ( : ) ==1 ) )  ,  A)  ; 

if  numel ( find (valids (:) ==1 )) >0 
for  i=l : length (valids ) 
if  valids (i)==l 

valid ( j ) =best_value ( i ) ; 
gain ( j ) =k  (i) ; 

p ( j ,  : )  =  [pi  (i)  p2(i)  p3(i)  p4(i)  p5(i)]; 
z ( j,  : )  =  [zl  (i)  z2(i)  z3(i)  z4(i)  z5(i)]; 

j=j+i; 

end 

end 

else 

aux2='F';  %  No  feasible  solution  exists. 

end 


case  9 

[best_value, mean_value, cl, k, zl, z2, z3, z4,pl,p2,p3,p4, c2]  = 
textread ( 'privl .txt ' , ' %f  %f  %c  %f  %f  %f  %f  %f  %f  %f  %f 
%s delimiter ' ,  '  '); 

[valids ] =t ext read ( ' validvector . txt ' , ' %u ' ) ; 


j=i; 

p=zeros (numel ( find (valids (:) ==1 ) ) , A) ; 
z  =  zeros (numel (find (valids ( : ) ==1 ) )  ,  A)  ; 

if  numel ( find (valids (:) ==1 )) >0 
for  i=l : length (valids ) 
if  valids (i)==l 

valid ( j ) =best_value ( i ) ; 
gain ( j ) =k  (i) ; 

p ( j  ,  : )  =  [pi  (i)  p2(i)  p3(i)  p4(i)]; 
z ( j,  : )  =  [zl  (i)  z2(i)  z3(i)  z4(i)]; 

j=j+i; 

end 

end 

else 

aux2='F';  %  No  feasible  solution  exists. 

end 


case  7 

[best_value, mean_value, cl , k, zl , z2 ,  z3 ,  pi ,  p2 ,  p3 ,  c2 ]  =  ... 
textread (' privl . txt ' , '%f  %f  %c  %f  %f  %f  %f  %f  %f  %n 
' delimiter ' ,  '  '); 

[valids ] =t ext read ( ' validvector . txt ' , ' %u ' ) ; 
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j=i; 

p=zeros (numel ( find (val ids {:) ==1 ) ) , A) ; 
z  =  zeros (numel (find (val ids ( : ) ==1 ) )  ,  A)  ; 

if  numel ( find (valids (:) ==1 )) >0 
for  i=l : length (valids ) 
if  valids (i)==l 

valid ( j ) =best_value ( i ) ; 
gain ( j ) =k  (i) ; 

p ( j ,  : )  =  [pi  (i)  p2(i)  p3(i)]; 
z ( j,  : )  =  [zl  (i)  z2  (i)  z3 (i)  ] ; 

j=j+i; 

end 

end 

else 

aux2='F';  %  No  feasible  solution  exists. 

end 
case  5 

[best_value, mean_value, cl , k, zl ,  z2 ,  pi ,  p2 ,  c2 ]  =  ... 

textread ( ' privl . txt ' ,  ' %f  %f  %c  %f  %f  %f  %f  %n  %s '  ,  .  .  . 
' delimiter ' ,  '  ' ) ; 

[valids ] =t ext read ( ' validvector . txt '  ,  ' %u ' ) ; 

j=i; 

p=zeros (numel ( find (valids (:) ==1 ) ) , A) ; 
z  =  zeros (numel (find (valids ( : ) ==1 ) )  ,  A)  ; 

if  numel ( find (valids (:) ==1 )) >0 
for  i=l : length (valids ) 
if  valids (i)==l 

valid ( j ) =best_value ( i ) ; 
gain ( j ) =k  (i) ; 
p ( j ,  : )  =  [pi  (i)  p2(i)]; 

z ( j,  :  )  =  [zl  (i)  z2  (i)  ]  ; 

j=j+i; 

end 

end 

else 

aux2='F';  %  No  feasible  solution  exists. 

end 


case  3 

[best_value, mean_value, cl , k, zl , pi , c2 ]  =  ... 

textread (' privl . txt %f  %f  %c  %f  %f  %n  %s delimiter ' ,  '  '); 
[valids ] =t ext read ( ' validvector . txt '  ,  ' %u ' ) ; 

j=i; 

p=zeros (numel ( find (valids (:) ==1 ) ) , A) ; 
z  =  zeros (numel (find (valids ( : ) ==1 ) )  ,  A)  ; 

if  numel ( find (valids (:) ==1 )) >0 
for  i=l : length (valids ) 
if  valids (i)==l 

valid ( j ) =best_value ( i ) ; 
gain ( j ) =k  (i) ; 
p ( j  ,  : )  =  [pi  (i)  ] ; 
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z(j,:)  =  [zl{i)]; 

j=j+i; 

end 

end 

else 

aux2='F';  %  No  feasible  solution  exists. 

end 


end 


m=A- 1 ; 

if  aux2=='T'  %  If  any  feasible  solution  exists. 

I=find (min (valid (:))) ;  %  Best  solution  of  each  iteration, 
if  auxx2==l 

zz=l;  %  Iterations  counter.  It  is  used  to  store  the  results. 


end 


%  Elements  to  store  the  results 
ac_gain= [ ] ; 

ac_p=zeros ( 9 , A) ; %  Number  of  possible 
%  of  poles=5 

ac_z  =  zeros ( 9 ,  A)  ;  %  Number  of  possible 
%  of  zeros=5 


controllers=9, max . number 
controllers=9, max . number 


ac_gain ( z  z ) =gain ( I ) ; 
for  1=1 : size (p, 2 ) 

ac_p ( z  z , 1 ) =p  ( 1 , 1 ) ; 
ac_z ( z  z , 1 ) =z (1,1) ; 

end 


% - Model  reduction - 

o, 

0 

%  Now  the  obtained  controller  is  reduced  to  a  m-th  order  model  by 
%  means  of  the  optimal  Hankel  norm  approximation. 
m=A- 1 ; 


if  m>=l 

sys=zpk (-ac_z (zz,  : ) , -ac_p (zz,  : ) , ac_gain (zz) )  ; 
sys=ss (sys) ; 
n=size (sys. A, 1) ; 
nu=size (sys, 2) ; 

[sysc,cinfo] =ncfmr (sys,  n)  ; 
sysch=hankelmr ( cinf o . GL,  m)  ; 
syschm=minreal ( inv ( sysch ( : , nu+1  . . . 

: end) ) *sysch ( : , 1 : nu) )  ; 


[ Z , P , K] =zpkdata ( syschm,  ' v ' ) ; 
z  z  =  z  z  +  1 ; 


%  Gain,  poles  and  zeros  of  the 
%  reduced  model. 


ac_gain ( z  z ) =K; 
for  1=1 : size (P ' , 2 ) 

ac_p  (zz,  1)  =-P  (1) 
ac_z (zz, 1) =-Z  (1)  '  ; 
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end 

z  z  =  z  z  +  1 ; 


end 


end 


if  m>=l 

clear  k  pi  p2  p3  p4  p5  zl  z2  z3  z4  z5  p  z  valids  valid  gain  . . . 

best_value  mean_value; 
delete  privl.txt 
delete  aagg.log 
delete  validvector.txt 

%  Update  the  CFG.txt  y  CFGCONTR.txt.  The  original  files  are  lost. 
File_Update; 
elseif  m<l 

auxx= 'F'; 

disp (' controllers  calculation  finished') 
for  h=l rnumel (find (ac_p ( : , 1) ~=0) ) 
ac_g= [ ] ; 
ac_poles= [ ] ; 
ac_zeros= [ ] ; 

for  1=1 rnumel (find (ac_p (h, : ) ~=0) ) 
ac_g (h) =ac_gain (h) ; 
ac_poles ( 1 ) =-ac_p (h, 1 )  ; 
ac_zeros (1) =-ac_z (h, 1) ; 

end 

G (h) =zpk (ac_zeros ( : ) , ac_poles ( : ) , ac_g (h) )  ; 
clear  ac_g  ac_poles  ac_zeros; 

end 

menuplot ; 

end 

auxx2=auxx2  +  l ; 

end 


%  B2AAGG  Function  to  generate  a  file  containing  the  bounds  variable,  in  an 

%  understable  format  for  the  optqft  program. 

o, 

0 

%  The  bounds  calculated  by  the  QFT  matlab  toolbox  has  the  following 

%  format: 

o, 

o 

%  -  each  column  corresponds  to  a  different  frequency 

%  -  if  there  are  n  analysed  points,  there  are  2n+2  rows 

%  -  by  default  there  are  73  phase  points  (  [ 0  : -5  : -3 60 ] ) 

%  -  the  first  n  values  correspond  to  the  minimum  limit 

%  (solid  line) 

%  -  the  next  n  values  correspond  to  the  maximum  limit 

%  (dashed  line) 

%  -  the  penultimate  value  corresponds  to  the  frequency 

%  -  the  last  value  indicates  the  type  of  bound 
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%  The  output  format  is: 

%  -  each  row  corresponds  to  a  different  frequency 

%  -  if  there  are  n  phase  analysed  points,  there  are  2n+2  columns 

%  -  the  first  value  corresponds  to  the  frequency 

%  -  the  next  value  corresponds  to  the  number  of  points  (73  by  default) 

%  -  the  next  n  values  correspond  to  the  minimum  limit 

%  -  the  next  n  values  correspond  to  the  maximum  limit 

%  Authors:  Mario  Garcia-Sanz,  Carlos  Molins 
%  1 . September . 2009 


function  []  =  b2aagg  (bounds.  Filename) 


if  (  ~isstr (Filename)  )  error  ('second  parameter:  file  name');  end 

[rows,  columns]  =  size  (bounds) ; 
numPhases  =  (rows-2)/2; 
b=bounds ( 1 : rows-2  ,:); 

b= [bounds ( rows-1 ,:) ;  numPhases*ones (1, columns)  ;  b]  ; 
b=b '  ; 

%  save  nombreFic  b  -ascii 

eval ( [ ' save  '  Filename  '  b  -ascii']); 


CHANGLE  Internal  utility  of  convexhl 

Authors:  Mario  Garcia-Sanz,  Carlos  Molins 
1 . September .2009 


function  ang=changle (pointl ,  point2) 

dx=real (point2) -real (pointl)  ; 
dy=imag (point2 ) -imag (pointl )  ; 
ang=atan2 (dy, dx) ; 
if(ang<0)  ang=ang+360;  end 


%  CHDIST  Internal  utility  of  convexhl.  Calculation  of  the  square  of  the 
%  distance. 

o, 

0 

%  Authors:  Mario  Garcia-Sanz,  Carlos  Molins 
%  1 . September . 2 00 9 


function  dist=chdist (pointl ,  point2) 

dx=real (point2) -real (pointl)  ; 
dy=imag (point2 ) -imag (pointl )  ; 
dist  =  dx*dx  +  dy*dy; 
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%  CHREMOVE  Internal  utility  of  convexhl .  It  removes  the  duplicated  points. 

o, 

0 

%  Authors:  Mario  Garcia-Sanz,  Carlos  Molins 
%  1 . September . 2 00 9 


function  v=chremove (complexvector ) 
v=  [  ]  ; 

[rows,  columns] =size (complexvector) ; 
for  i=l : rows 

[newrows,  columns] =size (v) ; 
repeated  =  0; 
for  j=l: newrows 

if  (equals (real (complexvector (i) ),  real (v { j )) )  &&  ... 

equals (imag (complexvector (1) ) , imag (v ( j ) ) )  ) 

repeated  =  1; 
break; 

end 

end 

if  (-repeated) 

v=[v;  complexvector ( 1 )]  ; 

end 

end 


%  CONVEXHL 

o, 

0 

o, 

o 

o, 

0 

o, 

o 

o, 

0 

o, 

o 

o, 

0 

o, 

o 

o, 

0 


Function  that  gives  the  convex-hull  of  a  set  of  complex  points 
,  but  calculated  in  the  db-phase  plane. 

complexvector:  Nxl  of  complex  points  vector. 

convex:  The  M  points  of  the  convex  polygon  that  contains  all 
the  points  arranged  anticlockwise. 

The  M  point  matches  with  the  1  in  order  to  close  the  polygon. 


%  Authors:  Mario  Garcia-Sanz,  Carlos  Molins 
%  1 . September . 2 00 9 


function  convex=convexhl (complexvector) 

[rows,  columns] =size (complexvector)  ; 

if  (columns  >  1)  error  ('it  has  to  be  a  Nxl  vector');  end 
if  (rows  <  3)  error  ('At  least  3  points  are  needed');  end 

vector=chremove (complexvector) ;  %  guitar  repetidos 

[rows,  columns] =size (vector) ; 

if  (rows  <  3)  error  ('At  least  3  different  points  are  needed');  end 

%  Mantain  two  complex  vectors: the  real  one  and  the  other  one  composed  of 
%  (phase  +  1  *  db) 
vectorReal  =  vector; 

[db, phase] =dbphase (vector) ; 
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vector  =  phase  +  i  *  db; 

%  The  calculations  are  made  in  "vector",  but  the  final  points  are  stored  in 
%  "vectorReal" 


%  Algorithm: 

o, 

0 

%  1.  Find  the  point  with  the  lowest  "y" 

%  2.  Exchange  between  the  first  element  and  the  lowest  one. 

%  The  first  position  is  i=l  point  of  the  convex-hull. 

%  3.  Find  the  i=i+l  of  the  convex-hull: 

%  measure  the  angle  between  the  i  point  and  the  rest  (with  respect  to 

%  the  X  axes) .  Take  into  account  that  the  convex-hull  angles  are 

%  always  increasing.  The  new  point  is  the  one  with  the  lowest  angle 

%  between  all  the  points  that  have  higher  angle  than  the  i  point. 


%  Calculate  the  minimum  y  value 

imin=l ; 

for  ii=2 : rows 

if (imag (vector (ii) ) <imag (vector (imin) ) )  imin=ii; 
elseif  ( imag (vector ( ii ) )  ==  imag (vector ( imin) ) ) 

%  Take  the  closest  one  to  the  left 

if (real (vector (ii) ) <real (vector (imin) ) )  imin=ii;  end 

end 

end 
%  swap 

tmp=vector ( 1 ) ;  vector ( 1 ) =vector ( imin) ;  vector (imin) =tmp; 

tmp=vectorReal (1) ;  vectorReal (1) =vectorReal (imin) ;  vectorReal (imin) =tmp; 

currentangle  =  0;  %  Angle  of  the  i-th  element,  the  i+l-th  element  must  be 

%  higher. 

%  For  each  element  i  of  the  convex-hull 
for  ii=l : rows 

%  Compare  with  the  remainder  looking  for  the  one  with  minimum  angle, 
if  (ii==l) 
jmin=2 ; 

angMin  =  changle (vector ( 1 ) ,  vector (2)); 
distMin  =  chdist (vector ( 1 )  ,  vector (2)); 
else 

jmin=l ; 

angMin  =  changle (vector ( ii ) ,  vector (1)); 
distMin  =  chdist (vector ( ii ) ,  vector (1)); 

end 

for  j=ii+l:rows 

ang  =  changle (vector ( ii ) ,  vector ( j) ); 
dist  =  chdist  (vector (ii),  vector (j)); 

%  Find  the  minimum  angle  between  all  that  beats  the  current, 
if  (ang  >=  currentangle) 
if (ang<angMin) 
jmin= j ; 
angMin  =  ang; 
distMin  =  dist; 
elseif  (ang  ==  angMin) 

%  Take  the  closest  one 
if  (dist<distMin) 
jmin= j ; 
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angMin  =  ang; 
distMin  =  dist; 

end 

end 

end 

end  %  for  j 

%  If  the  closest  one  is  1,  end 
if  (jmin==l) 

convex=vectorReal ( 1 : ii ) ; 

convex= [convex;  convex{l)];  %  close  the  polygon 
return 
else 

%  swap 

tmp=vector ( ii+1 ) ;  vector { ii+1 ) =vector { jmin) ;  vector ( jmin) =tmp; 
tmp=vectorReal (iitl) ;  vectorReal (ii  +  l) =vectorReal ( jmin)  ; 
vectorReal ( jmin) =tmp; 

end 

currentangle  =  angMin; 


end 


%  DBPHASE  Conversion  from  the  complex  plane  to  the  dB-phase  plane. 

o, 

0 

%  Authors:  Mario  Garcia-Sanz,  Carlos  Molins 
%  1 . September . 2009 


function  [db,  phase]  =  dbphase  (complexmatrix) 
db=20*logl0 (abs (complexmatrix) ) ; 

phase=l 8 0 *atan2 (imag (complexmatrix) , real (complexmatrix) ) /pi; 


o, 

0 

o, 

0 

o, 

0 

o, 

0 

o, 

0 

o, 

o 


EQUALS  Function  that  returns  1  if  vl  and  v2  are  approximately  equal. 


Authors:  Mario  Garcia-Sanz,  Carlos  Molins 
1 . September .2009 


function  isequal=equals (vl ,  v2  ) 
de it a= 0.0000001 ; 

isequal  =  (  (vl<v2+delta)  &  (vl>v2-delta)  ) ; 
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%  FILE_UPDATE  Function  to  update  the  CFG.txt  and  CFGCONTR.txt  files  during 
%  the  execution  of  the  Autoloopshaping_AG . 

o, 

0 

%  Authors:  Mario  Garcia-Sanz,  Carlos  Molins 
%  1 . September . 2 0 0 9 


function  File_Update ( ) 
fid  =  f open ( 'CFG.TXT r ')  ; 

font  =  f open ( ' outl . txt ' , ' wt ' ) ; 

a  =  textscan  ( f  id,  '  %d  %s  %s  %s%*[''\n]',  1); 
n  =  a{l}-2; 

fprintf ( f out ,  ' %d  %s  %s  %s\n ' , n, a { 2 } { 1 } , a { 3 } { 1 } , a { 4 } { 1 } ) ; 

for  i  =  1 : a { 1 } 

b  =  textscan ( fid,  ' %n  %s  %s  %s  %s  %s  %*[^\n]',  1); 
if  (i<=n) 

fprintf (font,  ' %d  %s  %s  %s  %s  %s  \n ' , b { 1 } , b { 2 } { 1 } , b { 3 } { 1 } , . . . 
b{4}  {l},b{5} {l},b{6} {!}) ; 

end 

end 

for  i  =  1 :  a { 1 } 

b  =  textscan ( fid,  ' %n  %s  %s  %s  %s 
if  (i<=n) 

fprintf  ( font ,  ' %d  %s  %s  %s  %s 

b{4} {l},b{5} {l},b{6} {!}) ; 

end 

end 

%  Parameters  Precision 
textLine  =  textscan ( fid,  ' %n  %s  %s  %s  %s  %s  %*[^\n]',  1); 
fprintf ( font ,  ' %d  %s  %s  %s  %s  %s  \n ', textLine { 1 }, textLine { 2 }{ 1 },..  . 

textLine{3} {!}, textLine { 4 } { 1 } , textLine { 5 } { 1 } , textLine { 6 } { 1 } )  ; 

%  Number  of  Individuals  in  the  population 
textLine  =  textscan  ( fid,  '  %d  %s  %*[''\n]',  1); 
fprintf ( font ,  ' %d  %s  \n ', textLine { 1 }, textLine { 2 }{ 1 }) ; 

%  Crossover  probability 

textLine  =  textscan  ( fid,  '  %n  %s  %s  %s  %*[''\n]',  1); 
fprintf ( font ,  ' %d  %s  %s  %s  \n ',  textLine { 1 },  textLine  {  2 }{ 1 },..  . 

textLine { 3 } { 1 } , textLine { 4 } { 1 } )  ; 

%  Mutation  probability 

textLine  =  textscan  ( fid,  '  %n  %s  %s  %s  %*["'\n]',  1); 
fprintf ( font ,  ' %d  %s  %s  %s  \n ', textLine { 1 }, textLine { 2 }{ 1 },..  . 

textLine { 3 } { 1 } , textLine { 4 } { 1 } )  ; 

%  Show  continuously  (1)  or  just  at  the  end  (0) 

textLine  =  textscan  ( fid,  '  %d  %s  %s  %s  %s  %s  %s  %s  %s  %*["'\n]',  1); 
fprintf ( font ,  ' %d  %s  %s  %s  %s  %s  %s  %s  %s  \n ',  textLine { 1 },  textLine { 2 }{ 1 },..  . 

textLine{3} {!}, textLine { 4 } { 1 } , textLine { 5 } { 1 } , textLine { 6 } { 1 } , . . . 
textLine{7}{l}, textLine { 8 } { 1 } ,  textLine { 9 }  { 1 } )  ; 


%s  %*  [-'Xn]  '  ,  1)  ; 

%s  \n' ,b{l},b{2} {l},b{3} {!}, . . . 
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%  Maximum  number  of  generations 

textLine  =  textscan ( f id,  ' %d  %s  %s  %s  %s  %*[^\n]',  1); 
fprintf ( f out ,  ' %d  %s  %s  %s  %s  \n textLine { 1 }, textLine { 2 }{ 1 },..  . 

textLine { 3 } { 1 } , textLine { 4 } { 1 } , textLine { 5 } { 1 } )  ; 

%  Minimum  desired  fitness 

textLine  =  textscan  ( fid,  '  %d  %s  %s  %s%*[''\n]',  1); 
fprintf ( font ,  ' %d  %s  %s  %s  \n ', textLine { 1 }, textLine { 2 }{ 1 },..  . 

textLine { 3 } { 1 } , textLine { 4 } { 1 } )  ; 

f close (font ) ; 
f close (fid) ; 

fid  =  f open ( 'CFGCONTR.TXT' ,  'r' )  ; 
font  =  f open ( ' out2 . txt ' , ' wt ' ) ; 

%Number  of  controller  fixed  poles 

a  =  textscan ( fid,  ' %d  %s  %s  %s  %s  %s  %s%*[^\n]',  1); 
fprintf ( font ,  ' %d  %s  %s  %s  %s  %s  %s\n ' , a { 1 } , a { 2 } { 1 } , a { 3 } { 1 } ,  . . . 

a{4}{l},a{5}{l},a{6}{l},a{7}{l}); 

%Number  of  parameters  that  correspond  to  the  controller  zeros 
a  =  textscan  ( fid,  '  %d  %s  %s  %s  %s  %s  %s  %s  %s  %s%*["'\n]',  1); 
n  =  a{l}-l; 

fprintf ( font ,' %d  %s  %s  %s  %s  %s  %s  %s  %s  %s  \n ' , n, a { 2 } { 1 } , a { 3 } { 1 } ,  .  .  . 
a{4}{l},a{5}{l},a{6}{l},a{7}{l},a{8}{l},a{9}{l},a{10}{l}); 

f close (font ) ; 
f close (fid) ; 

movefile  outl.txt  CFG. TXT 
movefile  out2.txt  CFGCONTR.TXT 

end 


%  FRECSN  Function  to  define  the  frequency  vector. 

o, 

0 

%  Authors:  Mario  Garcia-Sanz,  Carlos  Molins 
%  1 . September . 2009 

o, 

0 

function  w=frecsn() 

w=[0.01  0.05  0.1  0.2  1.0  5.0  10  20  30  40  50  100]; 


%  GRIDQFTN  Function  to  determine  the  number  of  points  of  each  plant 
%  parameter  (K,  Z,  P,  W,  D)  that  are  considered. 

o, 

o 

%  Authors:  Mario  Garcia-Sanz,  Carlos  Molins 
%  1 . September . 2009 
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function  nump=gridqf tn ( ) 
nump  =  3 ; 


o, 

0 


LIMITSN  Function  that  defines  the  range  of  variation  of  each  plant 
parameter . 

Plant  =  K  *  (1+s/z)  /  s  (l+s/p)  { l+2Ds/W+s^2/W^2 ) 

Authors:  Mario  Garcia-Sanz,  Carlos  Molins 
1 . September .2009 


function  [Kmin, Kmax, Zmin, Zmax, Pmin, Pmax, Wmin, Wmax, Dmin, Dmax] =limitsn ( ) 


Kmin  = 

0.2; 

Kmax  = 

2; 

Zmin  = 

0.5; 

Zmax  = 

0.75; 

Pmin  = 

1; 

Pmax  = 

10; 

Wmin  = 

5; 

Wmax  = 

6; 

Dmin  = 

CO 

o 

Dmax  = 

0.9; 

%  LOOPSHN  Plant  and  performance  specifications  definition  and  plant 
%  templates  and  QFT  bounds  creation.  Use  when  rectangular  or 

%  convex-hull  approximation  is  desired. 

o, 

o 

%  Authors:  Mario  Garcia-Sanz,  Carlos  Molins 
%  1 . September . 2009 


%  Nominal  Plant  Parameters  Values 
Knom  =  2 ; 

Znom  =  0.5; 

Pnom  =  10; 

Wnom  =  6; 

Dnom  =  0.8; 

%  Frequencies  to  consider 
w=f recsn; 

[m, numFrecs] =size (w) ; 

%  Templates  for  the  nominal  plant 

ptemplates=tmplnx (Knom,  Znom,  Pnom,  Wnom,  Dnom) ; 

%  Nominal  Plant 
nompt  =  1; 
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%%  Plant  templates  %% 
plottmpl (w, ptemplates , nompt )  ; 

%%  Bounds  calculation  %% 

Wl=[]  ; 

W2=[]  ; 

W3=  [  ]  ; 

W4=[]  ; 

W5=  [  ]  ; 

W6= [ ] ; 

W7=[]  ; 

W8=[]  ; 

W9=[]  ; 

R=0;  %  Delay 

%  Bound  type  1:  u-contour  in  6  dB  (robust  stability) 
wl  =  w; 

W1  =  10''(6/20);  %  20  log(x)  =  6 

boundl  =  sisobnds ( 1 , wl , Wl , ptemplates , R, nompt ) ; 

plotbnds (boundl) ; 

title ('Bound  type  1'); 

save  Wl 

%  Bound  type  7 :  Input  reference  tracking 
w7=[0.01  0.05  0.1  0.2  1.0  5.0  10.0]; 
numA  =  1 ; 

denA  =  conv  ([1  1],  conv  ([1  1],  [1/2  1])); 

inf  =  tf (numA, denA) ; 

numB  =  [1/0.35  1 ]  ; 

denB  =  conv  ([1/0.5  1],  [1/3  1 ] )  ; 

sup  =  tf (numB, denB) ; 

W7  =  [sup;  inf] ; 

boundl  =  sisobnds ( 7 , w7 , W7 , ptemplates , R, nompt ) ; 
plotbnds (boundl) ; 
title ('Bound  type  7'); 
save  W7 

%  Group  bounds 

bounds=grpbnds (boundl , boundl ) ; 

%  Bounds  intersection 
bounds=sectbnds (bounds)  ; 
plotbnds (bounds)  ; 
title  (' Bounds  Intersection'); 

%  'bnd_n'  file  generation  for  the  Genetic  Algorithm 
B2AAGG  (bounds,  'bnd_n') 

%  Frequency  array  for  the  loop-shaping 
wloop  =  logspace (-3, 3, 200) ; 
save  wloop 

%  Frequency  array  for  the  chksiso 
wcheck  =  logspace (-3, 3, 300) ; 
save  wcheck 

%  Nominal  Plant 

nump  =  Knom  *  [1/Znom  1]; 
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denp  =  conv  (  [10],  conv  ([1/Pnom  1],  [l/Wnom''2  2*Dnom/Wnom  1])  ); 

cp=tf (nump, denp)  ; 


Q, 

0 

O, 

0 

Q, 

0 


MENUPLOT  Menu  to  choose  the  task  to  do  with  the  results  obtained  from  the 
automatic  loop-shaping. 

Authors:  Mario  Garcia-Sanz,  Carlos  Molins 
1 . September .2009 


k=menu ( ' Choose  an  option ' ,  ' step ' ,  ' bode ' ,  ' nichols ' ,  ' nyquist '  ,  ' chksiso ' )  ; 

switch  k 

case  1 

for  h=l :numel (find (ac_p ( : , 1) ~=0) ) 
figure (2 ) 
step (cp*G (h) ) ; 
grid  on; 

set (gca,  ' FontSize ' ,  [5] ) 
axis  (  [0  70  0  200] ) 
hold  all 

end 
case  2 

for  h=l :numel (find (ac_p ( : , 1) ~=0) ) 
figure (2 ) 
bode (cp*G (h) , w) ; 
grid  on; 

set (gca,  ' FontSize ',  [5] ) 
axis([0.001  10000  -180  -45]) 
hold  all 

end 
case  3 

load  wloop; 

for  h=l :numel (find (ac_p ( : , 1) ~=0) ) 

Ipshape (wloop,  bounds,  cp,  G(h),  []); 

end 
case  4 

for  h=l :numel (find (ac_p ( : , 1) ~=0) ) 
figure (2 ) 
nyquist (cp*G (h) ) ; 
grid  on; 

set (gca,  ' FontSize  '  ,  [5] ) 
axis ( [-6  6  -3  3 ] ) 
hold  all 

end 
case  5 

load  wcheck; 
load  W1 
load  W7 

for  h=l :numel (find (ac_p ( : , 1) ~=0) ) 
if  ~isempty(Wl) 

chksiso ( 1 , wcheck, Wl,cp,0,G(h),l,l); 

title ([' bound  type  1  with  the  ' , num2str (h) , ' th  controller']); 

grid  on; 

end 

if  ~isempty(W2) 

chksiso (2 , wcheck, W2,cp,0,G(h),l,l); 

title ([' bound  type  2  with  the  ' , num2str (h) , ' th  controller']); 
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grid  on; 
end 

if  -isempty (W3 ) 

chksiso ( 3 , wcheck, W3,cp,0,G{h),l,l); 

title ([' bound  type  3  with  the  ' , num2str (h) , ' th  controller']); 

grid  on; 

end 

if  ~isempty (W4 ) 

chksiso ( 4 , wcheck, W4,cp,0,G{h),l,l); 

title ([' bound  type  4  with  the  ' , num2str (h) , ' th  controller']); 

grid  on; 

end 

if  ~isempty (W5 ) 

chksiso ( 5 , wcheck, W5,cp,0,G{h),l,l); 

title ([' bound  type  5  with  the  ' , num2str (h) , ' th  controller']); 

grid  on; 

end 

if  ~isempty (W6 ) 

chksiso ( 6 , wcheck, W6,cp,0,G{h),l,l); 

title ([' bound  type  6  with  the  ' , num2str (h) , ' th  controller']); 

grid  on; 

end 

if  ~isempty (W7 ) 

chksiso ( 7 , wcheck, W7,cp,0,G{h),l,l); 

title ([' bound  type  7  with  the  ' , num2str (h) , ' th  controller']); 

grid  on; 

end 

if  ~isempty (W8 ) 

chksiso ( 8 , wcheck, W8,cp,0,G{h),l,l); 

title ([' bound  type  8  with  the  ' , num2str (h) , ' th  controller']); 

grid  on; 

end 

if  ~isempty (W9 ) 

chksiso ( 9 , wcheck, W9,cp,0,G{h),l,l); 

title ([' bound  type  9  with  the  ' , num2str (h) , ' th  controller']); 

grid  on; 

end 

end 

end 


%  PLANTN  Function  to  calculate  the  plant  equation. 

o, 

0 

%  Planta  =  K  *  (1  +  s/z)  /  s  (l  +  s/p)  {l  +  2Ds/W+s"'2/W"'2) 

o, 

0 

%  Authors:  Mario  Garcia-Sanz,  Carlos  Molins 
%  1 . September . 2009 


function  g=plantn (K,  Z,  P,  W,  D,  w) 

%  w  is  a  row  vector 

g  =  K  .*  (l+j*w/Z)  ./  (  (j*w)  .*  {l+j*w/P)  .*  (1+2*0* j*w/W-w. *w/W^2)  ); 
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%  RECTANG  Function  that  makes  the  rectangular  approximation  of  the 

%  templates (the  smaller  rectangle  that  contains  the  original 

%  template) . 

g, 

0 

%  Each  column  is  a  templateand  each  row  indicates  a  new  point. 

%  NumPointSide  indicates  the  number  of  points  that  will  appear  in 

%  side  of  the  rectangle. 

o, 

0 

%  Authors:  Mario  Garcia-Sanz,  Carlos  Molins 
%  1 . September . 2 0 0 9 


function  tmplRect=rectang (templates, numPointSide) 

[numVertexes ,  numTemplates ] =size (templates)  ; 

[db, phase] =dbphase (templates) ; 

for  c=l : numTemplates 


mindb=l 0  0  0  0  0; 

o, 

0 

valor  minimo  en 

dB 

down=0 ; 

o, 

0 

indice  al  punto 

que  tiene 

mindb 

downphase=0 ; 

o, 

0 

fase  del  punto 

"aba jo" 

maxdb=- 100000; 

o, 

0 

valor  maximo  en 

dB 

up=0  ; 

o, 

0 

indice  al  punto 

que  tiene 

maxdb 

upphase=0 ; 

o, 

0 

fase  del  punto 

"arriba" 

minphase=l 0  0  0  0; 

o, 

0 

valor  minimo  en 

fase 

lef t=0 ; 

o, 

0 

indice  al  punto 

que  tiene 

minf ase 

dblef t=0 ; 

o, 

0 

dB  del  punto  "izq" 

maxphase=-l 0  0  0  0 

•  9- 

r  0 

valor  maximo  en 

fase 

right=0 ; 

o, 

0 

indice  al  punto 

que  tiene 

maxf ase 

dbright=0 ; 

o, 

0 

dB  del  punto  "der" 

for  f=l : numVertexes 

[db,  phase] 

= 

dbphase (templates (f, c) ) ; 

if  (db  >  maxdb)  maxdb  =  db;  up  =  f;  upphase  =  phase;  end 

if  (db  <  mindb)  mindb  =  db;  down  =  f;  downphase  =  phase;  end 

if  (phase  >  maxphase)  maxphase  =  phase;  right  =  f;  dbright  =  db;  end 

if  (phase  <  minphase)  minphase  =  phase;  left  =  f;  dbleft  =  db;  end 

end 

%  "right"  point  downwards  displacement  till  it  has  the  same  db  value  as 
%  the  "down"  point. 
displInDb  =  mindb  -  dbright; 

bottomrightpoint  =  templates (right, c)  *  10^ (displInDb/20) ; 

%  "up"  point  displacement  to  the  right  till  it  has  the  same  phase  value 
%  as  the  "right"  point. 

displInPhase  =  (maxphase  -  upphase) *pi/180; 
toprightpoint  =  templates (up, c)  *  exp ( j *displInPhase ) ; 

%  "left"  point  upwards  displacement  till  it  has  the  same  db  value  as 
%  the  "up"  point. 
displInDb  =  maxdb  -  dbleft; 

topleftpoint  =  templates  ( left ,  c)  *  10'' (displInDb/20)  ; 

%  "down"  point  displacement  to  the  left  till  it  has  the  same  phase  value 
%  as  the  "left"  point. 

displInPhase  =  (minphase  -  downphase) *pi/180; 
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bottomlef tpoint  =  templates (down, c)  *  exp { j *displInPhase ) ; 

%  Displacement  of  every  final  point  with  respect  to  the  previous  one. 

if  (numPointSide  <  2)  numPointSide  =  2;  end 

deltaInDb  =  (maxdb-mindb) / (numPointSide-l)  ; 

deltaInPhase  =  (maxphase-minphase) / (numPointSide-l )  ; 

deltaInPhase  =  deltaInPhase*pi/180; 

%  Finally,  all  the  points, 
for  f=l : numPointSide 

dispD=deltaInDb* (f-1) ; 
dispPh=deltaInPhase* (f-1)  ; 

tmplRect (f , c) =bottomrightpoint  *  10^{dispD/20); 
tmplRect (numPointSide+f , c) =toprightpoint  *  exp {- j*dispPh) ; 
tmplRect (2*numPointSide+f , c) =topleftpoint  *  10^ (-dispD/20) ; 
tmplRect (3*numPointSide+f , c) =bottomleftpoint  *  exp ( j*dispPh) ; 

end 


end 


%  TMPL_BNDS_AG  Plant  and  performance  specifications  definition  and  plant 
%  templates  and  QFT  bounds  creation.  Use  when  no  templates 

%  approximation  is  desired. 

o, 

0 

%  Authors:  Mario  Garcia-Sanz,  Carlos  Molins 
%  1 . September . 2009 


In  line  32  of  menuplot.m  write  " cp ( 1 , 1 , nompt ) "  where  "cp".If  not  it  will 
not  work 

%  Plant  definition 
c=l; 

for  K  =  linspace ( 0 . 2 , 2 , 3 ) 

for  Z  =  linspace ( 0 . 5 , 0 . 75 , 3 ) 
for  P  =  linspace ( 1 , 1 0 , 3 ) 

for  Wn  =  linspace ( 5 , 6 , 3 ) 

for  D  =  linspace ( 0 . 8 , 0 . 9 , 3 ) 
np(c, :)  =  K*[l/Z  1]; 

dp (c,  : )  =  conv { [ 1  0 ]  ,  conv { [ 1 /P  1 ]  ,  [ 1/Wn^2  2 *D/Wn  1 ] )  ) ; 

cp  (1,  1,  c)  =tf  (np  (c,  :  )  ,  dp  (c,  :  )  )  ; 
c=c+l ; 

end 

end 

end 

end 

end 

clear  K  Z  P  Wn  D  c; 

%  Nominal  plant 
nompt  =  187; 

%  Frequencies  of  interest 
w=f recsn; 
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%%  Templates  calculation  %% 

plottmpl (w, cp, nompt ) ; 

%%  Bounds  calculation  %% 

Wl=[]  ; 

W2=[]  ; 

W3=  [  ]  ; 

W4=[]  ; 

W5=  [  ]  ; 

W6= [ ] ; 

W7=[]  ; 

W8=[]  ; 

W9=[]  ; 

R=0;  %  Delay 

%  Bound  type  1:  u-contour  in  6  dB  (robust  stability) 
wl  =  w; 

W1  =  10^(6/20);  %  20  log(x)  =  6 

boundl  =  sisobnds ( 1 , wl , Wl , cp, R, nompt ) ; 
plotbnds (boundl) ; 
title ('Bound  type  1'); 
save  Wl 

%  Bound  type  7 :  Input  reference  tracking 
w7=[0.01  0.05  0.1  0.2  1.0  5.0  10.0]; 
numA  =  1 ; 

denA  =  conv  ([1  1],  conv  ([1  1],  [1/2  1])); 

inf  =  tf (numA, denA) ; 

numB  =  [1/0.35  1  ]  ; 

denB  =  conv  ([1/0.5  1],  [1/3  1 ] )  ; 

sup  =  tf (numB, denB) ; 

W7  =  [sup;  inf] ; 

boundl  =  sisobnds ( 7  ,  w7 , W7 , cp, R, nompt ) ; 
plotbnds (boundl) ; 
title ('Bound  type  7'); 
save  W7 

%  Group  bounds 

bounds=grpbnds (boundl , boundl ) ; 

%  Bounds  Intersection 
bounds=sectbnds (bounds)  ; 
plotbnds (bounds)  ; 
title (' Bounds  Intersection'); 

%  'bnd_n'  file  generation  for  the  Genetic  Algorithm 
B2AAGG  (bounds,  'bnd_n') 

%  Frequency  array  for  the  loop-shaping 
wloop  =  logspace (-3, 3, 200) ; 
save  wloop 

%  Frequency  array  for  the  chksiso 
wcheck  =  logspace (-3,3,300); 
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save  wcheck 

%  TMPLN  Function  to  calculate  the  plant  templates. 

%  Plant  =  K  *  (1  +  s/z)  /  s  (l  +  s/p)  {l  +  2Ds/W+s''2/W^2) 

o, 

0 

%  w  is  a  row  vector. 

%  Note  that  the  template  for  w=10  has  a  discontinuity  in  phase  and 

%  is  badly  calculated.  Due  to  this,  the  points  are  displaced. 

o, 

o 

%  Authors:  Mario  Garcia-Sanz,  Carlos  Molins 
%  1 . September . 2 0 0 9 


function  templates=tmpln ( ) 

[Kmin, Kmax, Zmin, Zmax, Pmin, Pmax, Wmin, Wmax, Dmin, Dmax] =limitsn; 
nump=gr idqf tn ; 
w=f recsn; 

complexmatrix  =  [ ] ; 

%  Generate  the  points 

for  K  =  vector (Kmin,  Kmax  ,nump) 

for  Z  =  vector (Zmin,  Zmax,  nump) 

for  P  =  vector (Pmin,  Pmax,  nump) 

for  W  =  vector (Wmin,  Wmax,  nump) 

for  D  =  vector (Dmin,  Dmax,  nump) 
line  =  plantn  (K,  Z,  P,  W, D, w) ; 
complexmatrix  =  [complexmatrix;  line]; 

end 

end 

end 

end 

end 

[mudo,  numTemplates] =size (w) ; 

%  Displace  w=10 

complexmatrix ( : , find (w==10) ) =complexmatrix ( : , find (w==10)  )  . *  exp  (- j*pi*3/2) ; 

templates= [ ] ; 

for  i=l : numTemplates 

convex=convexhl (complexmatrix ( : , i) ) ; 

%  There  may  be  dimension  problems  if  the  number  of  point  of  each 
%  convex-hull  does  not  match. 

[numVertexesAnt ,  mudo] =size (templates) ; 

[numVertexes ,  mudo] =size (convex) ; 
if  (numVertexesAnt>0 ) 

dif  =  numVertexes  -  numVertexesAnt; 
if  (dif  >  0) 

%  Complete  "templates" 

for  i=l:dif  templates= [templates ;  templates (numVertexesAnt ,:)] ;  end 
elseif  (dif  <  0) 

%  Complete  "convex" 

for  i=l:-dif  convex= [convex;  convex (numVertexes ,:)] ;  end 

end 

end 

templates  =  [templates  convex] ; 
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end 

%  TMPLNX  Function  to  calculate  the  Thompson-Nwokah  plant  templates,  adding 
%  the  nominal  plant  that  is  used  in  QFT .  This  function  is  similar 

%  to  TMPLN  with  the  difference  that  the  (Knom, Znom, Pnom, Wnom, Dnom) 

%  plant  is  added  as  the  first  of  the  templates. 

g, 

0 

%  *  TMPLN  is  used  to  calculate  the  area 

%  *  TMPLNX  is  used  as  template  in  QFT  (nomp  =  1) 

g, 

0 

%  w  is  a  row  vector. 

g, 

0 

g, 

0 

%  Authors:  Mario  Garcia-Sanz,  Carlos  Molins 
%  1 . September . 2 0 0 9 


function  tmpl=tmplnx (Knom,  Znom,  Pnom,  Wnom,  Dnom) 

%  The  first  plant  of  the  template 
w=f recsn; 

tmpl  =  plantn  (  Knom,  Znom,  Pnom,  Wnom,  Dnom,  w  ) ; 

%  Displacement  of  w=10 

tmpl (:, find (w==l 0 )) =tmpl (:, find {w==l 0 ) )  .*  exp  {-j*pi*3/2); 

%  The  remainder 
tmp=tmpln; 

%  Thompson-Nwokah  uses  a  rectangular  approximation  of  the  templates. 
tmp=rectang (tmp,  10); 

%  Remove  the  line  37  if  no  rectangular  approximation  is  desired. 

%  The  total 

tmpl  =  [tmpl;  tmp] ; 


%  VECTOR  Function  that  gives  vector  with  the  data  logarithmically  spaced 
%  out.  vmin  is  the  initial  value,  vmax  is  the  final  value  and  nump 

%  is  the  number  of  values  of  the  vector. 

g, 

o 

%  Authors:  Mario  Garcia-Sanz,  Carlos  Molins 
%  1 . September . 2009 

g, 

0 


function  v=vector (vmin, vmax, nump) 

V  =  logspace (loglO (vmin) , loglO (vmax) , nump); 


8.2.  C++  FILES 
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//  AAGG.h  Classes  definition  for  AAGG . cpp . 

// 

//  Authors:  Mario  Garcia-Sanz,  Carlos  Molins 
//  1 . September . 20 0 9 

// - 

#if  !  defined  _ ^AAGG_H _ 

# define  _ AAGG_H _ 

#include  <vector> 

#include  <fstream> 


//  Class  that  has  the  structure  of  each  solution  (genotype) 

//  A  Chromosome  is  a  set  of  booleans  (genes) 

//  (The  genes  are  grouped  to  form  parameters) 

// 

//  Index  of  the  first  gene  =  0 

// 

class  Chromosome 

{ 

private : 

bool  *  genes; 
int  numTotalGenes ; 

//  changes  one  gene 
void  change  (int  indGen) ; 

public : 

Chromosome  (  int  numG  ) ; 

Chromosome  (  const  Chromosome  &  father,  const  Chromosome  &  mother, 
int  Crosspoint  ) ; 

//  Crosspoint  =  0  — >  first  gene  of  the  father,  the  rest  of  the 
/ /  mother 

//  Crosspoint  =  numTotalGenes-1  — >  all  the  genes  of  the  father 
-Chromosome  (){  delete  []  genes;  } 

Chromosome  (  const  Chromosome  &  ) ; 

Chromosome  &  operator  =  (  const  Chromosome  &  ) ; 

//  gives  the  integer  value  represented  by  a  set  of  genes, 

//  supposing  that  the  one  with  the  lowest  index  is  the  less 
//  significative 

int  tovalue  (int  indGenInitial ,  int  numGenes)  const; 

//  mutates  genes  with  a  certain  probability 
void  mutate  ( ) ; 


//  Represents  a 
class  Parameter 
private : 


set  of  genes  that  forms  a  unity 

{ 
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int  numGenes; 

double  valueMin;  //  minimum  parameter  value  permitted 
double  valueMax;  //  maximum  parameter  value  permitted 
int  valueDecimalMax;  //  2**numGenes  -  1 

public : 

Parameter  (){}  //  necessary  to  use  later  std::vector 
Parameter  (  int  numG,  double  vMin,  double  vMax  ) ; 

//  By  default: 

//  Parameter  (  const  Parameter  &  )  ; 

//  Parameter  &  operator  =  (  const  Parameter  &  ) ; 

int  getNumGenes ( )  const  {  return  numGenes;  } 

//  given  a  set  of  genes,  it  translates  to  a  parameter  range 
double  tovalue  (  int  valueDecimal  )  const; 

//  necessary  to  use  later  std::vector 

bool  operator<  (  const  Parameter  &  p  )  const  {  return  true;  } 
bool  operator==  (  const  Parameter  &  p  )  const  {  return  false;  } 


//  Set  of  parameters,  shared  by  all  the  individuals: 
class  Parameters  { 


static  std::vector  <Parameter>  Parameterslist ; 
public : 

static  void  addParameter  (  int  numGenes,  double  valueMin, 
double  valueMax  ) ; 

static  const  Parameter  &  getParameter  (int  1) { 
return  Parameterslist [ 1 ]  ;  } 


static  int  getNumParameters  (){ 

return  Parameterslist . size ()  ;  } 

static  int  getNumTotalGenes  (); 


//  Phenotype:  coded  values  by  the  genotype  (saved  in  the  genes 
//  structure) :  one  value  for  each  parameter 

class  Phenotype  { 

private : 

double  *  values; 
public : 

Phenotype  (  const  Chromosome  &  ) ; 
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-Phenotype  (){  delete  []  values;  } 

Phenotype  (  const  Phenotype  &  ) ; 

Phenotype  &  operator  =  (  const  Phenotype  &  ) ; 

//  index  =  0  to  obtain  the  value  of  the  first  parameter 
double  getParameterValue  (int  index)  const  { 
return  values  [index];  } 


//  Utility  class  to  implement  the  scale  change  of  the  fitness 
//  function 
// 

//  The  fitness  function  f  has  two  problems: 

//  *  At  the  beginning  the  mediocre  individuals  disappear, 

//  because  they  are  clearly  worse  than  others.  Useful 

//  genetic  information  is  lost  (what  is  called  the 

//  "premature  convergence") . 

//  *  At  the  end  almos  all  the  individuals  are  similar. 

//  The  really  best  hardly  stand  out  ("random  walk"  effect) 

// 

//  The  solution  is  a  new  fitness  function  F=a*f+b,  that 

//  -  An  average  individual  of  f  has,  on  average,  one  offspring. 

//  F  must  mantain  this.  That  is  to  say,  f average==Faverage 

//  -  We  imposse  that  the  best  individual  of  F  has,  on  average, 

//  K  offsprings,  i.e.,  Fmax==K*Faverage 

// 

class  Scale  { 
private : 

static  const  int  K; 
static  double  a; 
static  double  b; 

public : 

//  calculates  the  new  a,b  given  certain  data  of  the 
//current  population 

static  void  recalculateScale  (  double  fitnessMin, 
double  fitnessMax,  double  fitnessAvg  ) ; 

//  evaluate  the  new  fitness  F 

static  double  getFitnessScaled  (  double  fitness  ) ; 

}; 


//  Represents  a  possible  solution  to  the  problem 
class  Individual  { 
private : 

Chromosome  *  genotype; 

Phenotype  *  phenotype; 

double  fitnessReal;  //  the  real 

mutable  double  f itnessofWork;  //  F=a*f+b 

//  the  f itnessofWork  is  unknown  at  the  birth  of  the 
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/ /  individual . 

//  It  os  known  when  all  the  generation  is  formed, 
void  evaluate  ();  //  calculate  the  real  fitness 

public : 

Individual  ( ) ; 

Individual  (  const  Individual  &  father, 

const  Individual  &  mother,  int  Crosspoint  ) ; 

//  Crosspoint  =  0  — >  first  gene  of  the  father,  the 
/ /  rest  of  the  mother 

//  Crosspoint  =  numTotalGenes-1  — >  all  the  genes  of 
//  the  father 
~Individual  ( )  ; 

Individual  (  const  Individual  &  ) ; 

Individual  &  operator  =  (  const  Individual  &  ) ; 

void  mutate  ( )  ; 

double  getFitnessReal  ()  const  {  return  f itnessReal ;  } 

double  getFitnessofWork  ()  const  { 
return  f itnessofWork;  } 

void  updateFitnessofWork  ()  const  {  f itnessofWork 

=  Scale :: getFitnessScaled  {  fitnessReal  )  ;  } 

bool  isValid  ()  const;  //  violates  or  not  the  possible 

//  constraints 


}; 


//  to  put  in  order  from  lower  to  higher,  with  the  best 
//  at  the  beginning 

bool  operator<  (  const  Individual  &  ind  )  const  {  return 
getFitnessReal ( )  >  ind . getFitnessReal {) ;  } 

bool  operator==  (  const  Individual  &  ind  )  const  {  return 
getFitnessReal ( )  ==  ind . getFitnessReal {) ;  } 

friend  std::ostream  &  operator<<  ( std : : ostream  &  os,  const 
Individual  &); 


std::ostream  &  operator<<  ( std :: ostream  &  os,  const  Individual 

&  indiv) ; 


//  Utility  class  to  implement  the  selection  by  means  of  the 
//  "stochastic  remainder" 

class  SRemainder  { 

private : 

//  each  element  of  the  array  is  an  index  between  0  and 
/ /numIndividuals-1 
int  *  arrayindexes ; 
int  numindividuals ; 

int  unSelected;  //  number  of  individuals  that  are  to 
//  be  selected 

SRemainder  (  const  SRemainder  &  ) ; 

SRemainder  &  operator  =  (  const  SRemainder  &  ) ; 
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public : 

SRemainder  (  const  std::vector  <Individual>  & 

Individualslist ,  double  avgFitnessofWork  ); 
~SRemainder  (); 

int  get Index  (); 

}; 


//  Represents  the  set  of  individuals  put  in  order  by  their 
//  fitness  value  (from  best  to  worst) 


class  Population  { 
private : 

std::vector  <Individual>  Individualslist; 
double  sumFitnessReal ;  //  sum  fitnessReal  of 

/ /  each  individual 

double  maxFitnessReal ;  //  maximum  fitnessReal  value 

//  that  is  in  the  population 
double  minFitnessReal ;  //  minimum  fitnessReal  value 

//  that  is  in  the  population 
double  sumFitnessofWork;  //  sum  f itnessofWork  of 

/ /  each  individual 

mutable  SRemainder  *  sremainder; 


Population  (  const  Population  &  ) ; 

Population  &  operator  =  (  const  Population  &  ) ; 


public : 

Population  ();  //  empty  population 

Population  (  int  numindividuals  );//  random  population 
-Population  (); 


void  addindividual  (  const  Individual  &  ) ; 
//  calculates  the  f itnessofWork 
void  endof aGeneration  (); 


}; 


const  Individual  &  select  ()  const; 

double  getFitnessRealAverage  ()  const  {  return 

sumFitnessReal  /  (double)  Individualslist . size () ;  } 

double  getMaxFitnessReal  ()  const  {  return 
maxFitnessReal;  } 

double  getMinFitnessReal  ()  const  {  return 
minFitnessReal;  } 

double  getBestFitnessReal  ()  const  {  return 

Individualslist [ 0 ] . getFitnessReal ( ) ;  } 

double  getFitnessofWorkAverage  ()  const  {  return 

sumFitnessofWork  /  (double)  Individualslist . size () ;  } 

const  Individual  &  getBestIndividual  ()  const  {  return 
Individualslist  [  0 ] ;  } 


std::ostream  &  operator<<  ( std : : ostream  &  os,  const 

Population  &  pob) ; 
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//  Utility  class  that  represents  the  Genetic  Algorithm  itself 
class  GeneticAlgorithm  { 
private : 

Population  *  thePopulation; 

Population  *  theNewPopulation;  //  next  generation 

int  numGenerations ; 

bool  showContinuous ; 

int  numGenerationsMax; 

double  minFitnessRealDesired; 

GeneticAlgorithm  (  const  GeneticAlgorithm  &  ) ; 
GeneticAlgorithm  &  operator  =  (  const  GeneticAlgorithm  &  ) ; 

const  Individual  &  select  (); 

void  cross  (const  Individual  &  father,  const 
Individual  &  mother) ; 
void  mutate  (Individual  &  indiv) ; 
void  passGeneration  ( )  ; 

bool  stop  0;  //  indicates  the  end  of  the  algorithm, 

void  show  (  std::ostream  &  os  )  ; 

public : 

GeneticAlgorithm  ( ) ; 

~GeneticAlgorithm  ( )  ; 

void  initialize  (  int  numindividuals ,  double  pCross,  double 
pmutate,  bool  showContinuous,  int  maxGenerations , 
double  minFitnessReal  ) ; 
void  mainAlgorithm  (  std::ostream  &  os  )  ; 


}; 

#endif 


Basic  Genetic  Algorithm . C++  adaptation  of 
Goldberg  basic  code  [23] . 


/ /AAGG . cpp 

// 

// 

// 

//  Authors:  Mario  Garcia-Sanz,  Carlos  Molins 
//  1 . September . 20 0 9 

// - 


the 


#include  <math.h>  //  for 

#include  <time.h>  //  for 

#include  <algorithm>  //  for 

#include  <float.h>  //  for 

#include  "random. h" 

#include  "aagg.h" 

/ /  Problem  dependent  parameters 
int  numIndividualsInPopulation; 
double  pMutate; 
double  pCross; 


floor  ( ) 
time ( ) 
sort  ( ) 
DBL_MAX 
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//  External  evaluation  function 

extern  double  externalEvaluation  (  const  Phenotype  &  )  ; 

//  Function  to  check  the  external  validity 
//  (no  constraint  violation) 

extern  bool  isValidExternalSolution  (  const  Phenotype  &  ) ; 


Chromosome :: Chromosome  (  int  numG  ){ 

genes  =  new  bool  [numTotalGenes  =  numG  ] ; 

//  random  initialization 

for  (  int  i  =  0;  i  <  numG;  i++  )  genes  [i]  =  flip  (0.5); 


Chromosome :: Chromosome  (  const  Chromosome  &  father,  const  Chromosome  & 

mother,  int  Crosspoint  ) { 

//  Crosspoint  =  0  — >  father's  first  gene,  the  rest  of  the  mother 
//  Crosspoint  =  numTotalGenes-1  — >  all  the  gens  of  the  father 
genes  =  new  bool  [  numTotalGenes  =  father . numTotalGenes  ] ; 
for  (  int  i  =  0;  i  <=  Crosspoint;  i++  )  genes  [i]  =  f ather . genes [ i ] ; 
for  (  int  i  =  Crosspoint  +1;  i  <  numTotalGenes;  i++  )  genes  [i]  = 
mother . genes [ i ] ; 

} 

Chromosome :: Chromosome  (  const  Chromosome  &  crom  ){ 
numTotalGenes  =  crom. numTotalGenes; 
genes  =  new  bool  [numTotalGenes  ] ; 

for  (int  i  =  0;  i  <  numTotalGenes;  i++)  genes [i]  =  crom . genes [i]; 


Chromosome  &  Chromosome :: operator  =  (  const  Chromosome  &  crom  ) { 
if  (  &crom  !=  this  ) { 
delete  []  genes; 

numTotalGenes  =  crom. numTotalGenes; 
genes  =  new  bool  [numTotalGenes  ] ; 

for  (int  i  =  0;  i  <  numTotalGenes;  i++)  genes [i]  =  crom . genes [ i ] ; 

} 

return  *this; 


int  Chromosome :: tovalue  (int  indGenInitial ,  int  numGenes)  const  [ 
int  indGenFinal  =  indGenInitial  +  numGenes  -  1; 
if  (indGenFinal  >=  numTotalGenes)  throw; 

int  value  =  0  ; 
int  powerofC  =  1; 

for  (int  i=indGenInitial ;  i  <=  indGenFinal;  i++) { 
if  (genes  [i])  value  +=  powerof2; 
powerofC  *=  2; 

} 

return  value; 


void  Chromosome: :mutate  () { 

for  (int  i=0;  i  <  numTotalGenes;  i++) { 
if  (  flip  (pMutate  )  )  change  (i) ; 

} 

} 
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void  Chromosome :: change  (int  indGen) { 
genes [indGen]  =  !  genes [indGen] ; 

} 


Parameter :: Parameter  (  int  numG,  double  vMin,  double  vMax  ) 

: numGenes (numG) ,  valueMin (vMin)  ,  valueMax (vMax)  { 
int  scale  =  1; 

scale  =  scale  <<  numGenes;  //  2  to  the  power  of  numGenes 
valueDecimalMax  =  scale  -  1; 


double  Parameter :: tovalue  (int  valueDecimal )  const  { 

//  translates  from  the  range  [0,  2**L  -  1]  to  the  range  [min,  max] 
double  range  =  valueMax  -  valueMin; 

return  valueMin  +  ((double)  valueDecimal)  *  range  /  valueDecimalMax; 


//  Parameters  list 

std::vector  <Parameter>  Parameters :: Parameterslist ; 

void  Parameters :: addParameter  (  int  numGenes,  double  valueMin,  double 

valueMax  ) { 

Parameter  param  (numGenes,  valueMin,  valueMax); 

Parameterslist . push_back  (param) ; 


int  Parameters :: getNumTotalGenes  (){ 
int  numGenes  =  0; 

for  (int  1=0;  1  <  getNumParameters ( ) ;  i++) { 

numGenes  +=  getParameter ( 1 ) . getNumGenes ( ) ; 

} 

return  numGenes; 

} 

Phenotype :: Phenotype  (  const  Chromosome  &  genotype  ){ 
int  numParameters  =  Parameters :: getNumParameters () ; 
values  =  new  double  [numParameters]; 
int  IndGenlnitial  =  0; 

for  (int  1=0;  1  <  numParameters;  i++) [ 

int  numGenes  =  Parameters :: getParameter ( 1 ) . getNumGenes () ; 
int  valueDecimal  =  genotype . tovalue  (IndGenlnitial,  numGenes); 
values[i]  =  Parameters :: getParameter ( 1 ). tovalue (valueDecimal )  ; 
IndGenlnitial  +=  numGenes; 

} 


Phenotype :: Phenotype  (  const  Phenotype  &  f  ){ 

int  numParameters  =  Parameters :: getNumParameters ()  ; 
values  =  new  double  [numParameters]; 

for  (int  1=0;  1  <  numParameters;  i++)  values [1]  =  f. values [1]; 

} 

Phenotype  &  Phenotype :: operator  =  (  const  Phenotype  &  f  ){ 
if  (  &f  ! =  this  ) [ 

int  numParameters  =  Parameters :: getNumParameters () ; 
delete  []  values; 

values  =  new  double  [numParameters]; 

for  (int  1=0;  1  <  numParameters;  i++)  values [1]  =  f. values [1]; 

} 

return  *this; 
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const  int  Scale: :K  =2;  //  it  must  be  >1 

double  Scale: :a  =  1; 
double  Scale: :b  =  0; 

void  Scale :: recalculateScale  (  double  fitnessMin,  double  fitnessMax, 

double  fitnessAvg  ) { 

if  (fitnessMin  >  (  (K*f itnessAvg-f itnessMax) / (K-l))  ) { 

double  delta  =  fitnessMax  -  fitnessAvg; 
a  =  (K-l)  *  fitnessAvg  /  delta; 

b  =  fitnessAvg  *  (fitnessMax  -  K  *  fitnessAvg)  /  delta; 

} 

else  { 

double  delta  =  fitnessAvg  -  fitnessMin; 
a  =  fitnessAvg  /  delta; 

b  =  -  fitnessMin  *  fitnessAvg  /  delta; 

} 

} 


double  Scale :: getFitnessScaled  (  double  fitness  ){ 
return  a  *  fitness  +  b; 

} 


Individual :: Individual  (){ 

//  create  a  random  individual 

int  numTotalGenes  =  Parameters :: getNumTotalGenes () ; 
genotype  =  new  Chromosome  (numTotalGenes); 
phenotype  =  new  Phenotype  (*genotype); 
evaluate  (); 

} 


Individual :: Individual  (  const  Individual  &  father, 

const  Individual  &  mother,  int  Crosspoint 

)  { 

//  Crosspoint  =  0  — >  first  gen  of  father,  the  rest  of  mother 
//  Crosspoint  =  numTotalGenes-1  — >  all  the  gens  of  father 
genotype  =  new  Chromosome  ( *father . genotype,  ^mother . genotype. 

Crosspoint) ; 

phenotype  =  new  Phenotype  (*genotype); 
evaluate  (); 

} 


Individual :: ~Individual  (){ 
delete  genotype; 
delete  phenotype; 

} 


Individual :: Individual  (  const  Individual  &  indiv  ){ 
genotype  =  new  Chromosome  (*indiv. genotype) ; 
phenotype  =  new  Phenotype  (*indiv. phenotype) ; 
fitnessReal  =  indiv . fitnessReal ; 
f itnessofWork  =  indiv . fitnessofWork; 

} 


Individual  &  Individual :: operator  =  (  const  Individual  &  indiv  ){ 
if  (  Sindiv  !=  this  ) { 

*genotype  =  *indiv. genotype; 
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*phenotype  =  *indiv. phenotype; 
fitnessReal  =  indiv . f itnessReal  ; 
f itnessofWork  =  indiv . fitnessofWork; 

} 

return  *this; 


void  Individual :: evaluate  (){ 

//  calculate  fitness 

fitnessReal  =  externalEvaluation  { *phenotype) ; 
fitnessofWork  =  0.0; 

} 

void  Individual : :mutate  (){ 
genotype->mutate ( )  ; 
delete  phenotype; 

phenotype  =  new  Phenotype  (*genotype); 
evaluate  (); 

} 


bool  Individual :: isValid  ()  const  { 

//  it  fulfills  or  not  the  constraints 

return  isValidExternalSolution  (  *phenotype  ) ; 

} 


std : : ostream  &  operator<<  ( std : : ostream  &  os,  const  Individual  &  indiv) { 
int  numParameters  =  Parameters :: getNumParameters {) ; 

os  «  "  ["; 

for  (int  i  =  0;  i  <  numParameters;  i++) { 

os  <<  indiv . phenotype->getParameterValue  (i)  <<  "  "; 

} 

return  os  <<  "]"  <<  std::endl; 


SRemainder : : SRemainder  (  const  std::vector  <Individual>  &  Individualslist , 

double  avgFitnessofWork  ) { 

//  based  on  "stochastic  remainder  without  replacement" 

numindividuals  =  Individualslist . size {) ; 
arrayindexes  =  new  int  [numindividuals]; 
unSelected  =  numindividuals; 

//  fill  the  array 

//  each  element  of  the  array  is  an  index  between  0  and  numIndividuals-1 

double  *  fractions  =  new  double  [numindividuals]; 
int  ind  =  0;  //  one  index  to  the  array  of  indexes 

//  assign  the  integer  number  of  copies 
for  (int  i  =  0;  i  <  numindividuals;  i++) { 

double  hopedOff springs  =  Individualslist [ i ]. getFitnessofWork ( )  / 

avgFitnessofWork; 

int  Fixedoff springs  =  (int)  floor  (hopedOff springs ) ; 

fractions [i]  =  hopedOff springs  -  (double)  Fixedoff springs ; 

for  (int  j  =  0;  j  <  Fixedoff springs ;  j++)  arrayindexes [ind++]  =  i; 

} 

//  assign  the  fractional  number  of  copies 
int  i  =  0; 

while  (ind  <  numindividuals) { 


102 


Advanced  Quantitative  Robust  Control  Engineering:  Ref.  BOARD  Award  #  FA8655-08-1-3027 

New  Solutions  for  Automatic  Loop-shaping  for  SISO  Ref.  Grant  #  083027 

and  MIMO  Systems.Part  1:  SISO  Systems.  Date:  l.September.2009  Page  103 


if  (  fractions  [i]  >  0  ) { 

bool  win  =  flip  (  fractions  [i]  ) ; 

if  (win) { 

arraylndexes [ ind++ ]  =  i; 

//  not  to  win  again 
fractions  [i]  =  0 ; 

} 

} 

i  +  +; 

if  ( i>=numlndividuals )  i=0; 

} 

delete  []  fractions; 


SRemainder: : ~SRemainder  () { 
delete  []  arraylndexes; 

} 

int  SRemainder :: getindex  (){ 

//  Give  all  the  array  indexes  one  by  one.  The  already  given 
//  indexes  are  at  the  top  and  the  rest  at  the  bottom. 

//  It  gives  a  random  index  from  the  array  of  indexes  between  the 
//  indexes  not  given  yet. 
if  (  unSelected  ==  0  )  throw; 

int  sel  =  getRandomIndex  (0,  unSelected  -  1); 
int  result  =  arraylndexes  [sel]; 

//  move  to  the  begining  the  last  index, 
arraylndexes [ sel ]  =  arraylndexes [unSelected  -  1]; 
unSelected — ; 
return  result; 

} 

Population :: Population  (){ 
sumFitnessReal  =  0.0; 
maxFitnessReal  =  0.0; 
minFitnessReal  =  DBL_MAX; 
sumFitnessofWork  =  0.0; 
sremainder  =  0; 


Population :: Population  (  int  numindividuals  ){ 
sumFitnessReal  =  0.0; 
maxFitnessReal  =  0.0; 
minFitnessReal  =  DBL_MAX; 

for  (int  i  =  0;  i  <  numindividuals;  i++) { 

Individual  ind;  / /  random 
Individualslist . push_back  (ind) ; 
double  fit  =  ind . getFitnessReal ( )  ; 
sumFitnessReal  +=  fit; 

if  (fit  >  maxFitnessReal)  maxFitnessReal  =  fit; 
if  (fit  <  minFitnessReal)  minFitnessReal  =  fit; 

} 

std::sort  (  Individualslist . begin ()  ,  Individualslist . end ( )  ); 

sremainder  =  0; 

endof aGeneration  ();  //  calculate  f itnessofWork 


Population :: -Population  (){ 

if  (sremainder)  delete  sremainder; 

} 
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void  Population :: addindividual  {  const  Individual  &  ind  ){ 
Individualslist . push_back  (ind) ; 
double  fit  =  ind . getFitnessReal ( )  ; 
sumFitnessReal  +=  fit; 

if  (fit  >  maxFitnessReal )  maxFitnessReal  =  fit; 
if  (fit  <  minFitnessReal )  minFitnessReal  =  fit; 
std::sort  (  Individualslist . begin ()  ,  Individualslist . end ( )  ); 

} 


void  Population :: endofaGeneration  (){ 

Scale :: recalculateScale  (minFitnessReal,  maxFitnessReal, 
getFitnessRealAverage ( )  ); 

typedef  std :: vector<Individual> :: iterator  iter; 
sumFitnessofWork  =  0.0; 

for  (iter  i  =  Individualslist . begin ()  ;  i  !  = 
Individualslist . end  0 ;  ++i  ){ 
i->updateFitnessofWork  ( )  ; 

sumFitnessofWork  +=  i->getFitnessofWork ( ) ; 

} 

} 


const  Individual  &  Population :: select  ()  const  { 

//  Stochastic  Remainder  method 

if  (sremainder  ==  0  )  sremainder  =  new  SRemainder  (  Individualslist, 
getFitnessofWorkAverage  ( )  ); 

return  Individualslist!  sremainder->getlndex ( )  ] ; 


std::ostream  &  operator<<  ( std : : ostream  &  os,  const  Population  &  pob) { 
if  (  !  pob . getBestIndividual  () .isValidO  )  os  << 

"this  is  not  a  feasible  solution"  <<  std::endl; 
os  <<  "Best  value  =  "  <<  pob . getBestFitnessReal  ()  <<  std::endl 
<<  "Mean  value  =  "  <<  pob . getFitnessRealAverage  ()  <<  std::endl 
<<  "Parameters  =  "  <<  pob . getBestIndividual  ()  <<  std::endl  << 
std : : endl ; 

//  To  obtain  the  evolution  graphics 

std :: of stream  ofs  ("privl.txt",  std : : ios : : app)  ; 

ofs  <<  pob . getBestFitnessReal  ()  <<  "\t"  <<  pob . getFitnessRealAverage  () 
<<  "\t"  <<  pob . getBestIndividual  ()  <<  std::endl; 
ofs . close  ( ) ; 

//  To  save  a  vector  with  1  for  feasible  solution  and  0  for  the 
//  non-feasible  ones 

std :: of stream  ovv  ("validvector.txt",  std : : ios : : app) ; 

if  (  !  pob . getBestIndividual  () .isValidO  )  ovv  <<  "0"  <<  std::endl; 

if  (  pob . getBestIndividual  () .isValidO  )  ovv  <<  "1"  <<  std::endl; 

ovv . close ( )  ; 

return  os; 


GeneticAlgorithm : : GeneticAlgorithm  () { 
thePopulation  =  0; 
theNewPopulation  =  0; 
numGenerations  =  0; 
showContinuous  =  false; 
numGenerationsMax  =  0; 
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minFitnessRealDesired  =  0 ; 

} 


GeneticAlgorithm : : ~GeneticAlgorithm  {) { 

if  (thePopulation)  delete  thePopulation; 
if  (theNewPopulation)  delete  theNewPopulation; 

} 


void  GeneticAlgorithm: : initialize  {  int  numindividuals ,  double  pCr,  double 
pMt, 

bool  showCt,  int  ngMax, 

double  mf  ) { 

if  (  Parameters :: getNumParameters  {)  ==  0  )  throw; 

//  even  number  of  individuals 
if  (  numindividuals  %  2  )  throw; 

//  global  variables 

numIndividualsInPopulation  =  numindividuals; 
pMutate  =  pMt; 
pCross  =  pCr; 

/ /  random  numbers 
srand  (  time  (NULL)  )  ; 

//  first  generation 

thePopulation  =  new  Population  (numindividuals); 
theNewPopulation  =  new  Population  (); 
numGenerations  =  1; 

//  member  data 
showContinuous  =  showCt; 
numGenerationsMax  =  ngMax; 
minFitnessRealDesired  =  mf; 

} 

void  GeneticAlgorithm: :mainAlgorithm  (  std::ostream  &  os  )  { 
while  (  !  stop ( )  ) { 

//  suppose  a  even  number  of  individuals 

for  (  int  i  =  0;  i  <  numIndividualsInPopulation;  i=i+2  ) { 
const  Individual  &  father  =  select (); 
const  Individual  &  mother  =  select (); 
cross  (  father,  mother  ) ;  //  mutation 

} 

if  (showContinuous)  show(os); 
passGeneration ( ) ; 

} 

if  (!  showContinuous)  show (os); 

} 


const  Individual  &  GeneticAlgorithm :: select  (){ 
return  thePopulation->select  (); 

} 


void  GeneticAlgorithm: : cross  (const  Individual  &  father, 

const  Individual  &  mother) { 
int  numGenes  =  Parameters ::getNumTotalGenes(); 


int  Crosspoint  =  numGenes  -  1; 
//  Crosspoint  =  0  — >  first  gen 
//  Crosspoint  =  numTotalGenes-1 
if  (  flip  (pCross)  )  Crosspoint 
/ /  offsprings 

Individual  offspringl  (  father. 
Individual  offspring2  (  mother. 


of  father,  the  rest  of  mother 
— >  all  the  gens  of  father 
=  getRandomIndex  (0,  numGenes  -  2); 

mother.  Crosspoint  ) ; 
father.  Crosspoint  ) ; 
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//  mutation 

of f springl .mutate ( ) ; 

of f spring2 .mutate  ( ) ; 

//  save 

theNewPopulation->addIndividual  {  offspringl  ) ; 
theNewPopulation->addIndividual  {  offspring2  ) ; 


void  GeneticAlgorithm: :mutate  (Individual  &  indiv) { 
indiv. mutate ( ) ; 

} 


void  GeneticAlgorithm: :passGeneration  {){ 
delete  thePopulation; 
thePopulation  =  theNewPopulation; 

thePopulation->endof aGeneration  { ) ;  //  calculate  f itnessofWork 

theNewPopulation  =  new  Population  (); 
numGene rat ions ++; 


bool  GeneticAlgorithm: : stop  (){ 

return  (numGenerations  >=  numGenerationsMax) 
(thePopulation->getBestFitnessReal  { ) 


>=  minFitnessRealDesired)  ; 


void  GeneticAlgorithm: : show  (  std::ostream  &  os  ) { 

os  <<  "Generation  "  <<  numGenerations  <<  std::endl  <<  *thePopulation; 

} 


//  BOUNDS. h  Classes  definition  to  represent  the  QFT  bounds  and  to 

//  calculate  the  cost  of  a  controller.  The  fitness  function 

//  is  based  on  Thompson  [39] . 

// 

//  Authors:  Mario  Garcia-Sanz,  Carlos  Molins 
//  1 . September . 20 0 9 

// - 

#if  [defined  _ BOUNDS_H _ 

# define  _ BOUND S_H _ 

#include  <complex> 

#include  <vector> 

#include  <algorithm>  //  for  sort() 

#include  <iostream> 

#include  "aagg.h" 


class  qcomplex  { 
private : 

std : : complex<double>  c; 
public : 

qcomplex  (double  re  =  0,  double  im  =  0)  :  c(re,im){} 

qcomplex  ( std : : complex<double>  cc)  :  c(cc){} 
-qcomplex  ( ) { } 
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}; 


//  qcomplex  (  const  qcomplex  &  ) ; 

//  qcomplex  &  operator  =  (  const  qcomplex  &  ) ; 

double  realO  const  {  return  c.realO;  } 
double  imagO  const  {  return  c.imagO;  } 
double  abs()  const  {  return  std::abs(c);  } 
double  argO  const  {  return  std::arg(c);  } 
qcomplex  conjO  const  {  return  std:  :  con  j  (c)  ;  } 

operator  std : : complex<double>  ()  {  return  c;  } 

operator  std : : complex<double>  ()  const  {  return  c;  } 

//  to  keep  arranged  the  zeros  of  a  polynomial 
//  from  lower  to  higher,  with  the  most  significative  at  the 
/ /  beginning . 

bool  operator<  (  const  qcomplex  &  c2  )  const  { 
return  real()  <  c2.real();  } 
bool  operator==  (  const  qcomplex  &  c2  )  const  { 

return  real()  ==  c2.real()  &&  imag { )  ==  c2.imag{);  } 
friend  std::istream  &  operator>>  ( std : : istream  &  is, 
qcomplex  &  qc) ; 


std::istream  &  operator>>  ( std :: istream  &  is,  qcomplex  &  qc) { 
return  is  >>  qc.c; 

} 

//  Represents  a  numerator  or  a  denominator  of  a  rational  transfer 
//  function  in  the  form  (s+zl) (s+z2) . . .  where  zl,  z2,  . . .  may  be 

//  complex. 

// 

//  The  list  is  kept  arranged  from  the  most  to  the  less  sinificative 
//  in  order  to  select  the  first  n  zeros. 

class  Polynom  { 

std::vector  <qcomplex>  Zeroslist; 
int  numOriginZeros ; 

public : 


Polynom  () {  numOriginZeros  =  0;  } 

//  By  default: 

//  Polynom  (  const  Polynom  &  ) ; 

//  Polynom  &  operator  =  (  const  Polynom  &  ) ; 

void  addZero  (  qcomplex  c  ) ; 

int  getNumZeros  ()  const  {  return  Zeroslist . size  () ;  } 

int  getNumZerosNotZero  ()  const  { 

return  Zeroslist . size ( )  -  numOriginZeros;  } 
qcomplex  getZero  (int  i)  const  {  return  Zeroslist [ i ] ;  } 

qcomplex  getZeroNotZero  (int  i)  const  { 

return  Zeroslist [ i+numOriginZeros ] ;  } 

qcomplex  evaluatein  (  double  frequency  )  const; 

friend  Polynom  operator  *  (  const  Polynom  &  pi, 

const  Polynom  &  p2  ) ; 
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Polynom  operator  *  (  const  Polynom  &  pi, 

const  Polynom  &  p2  ) ; 


//  Represents  a  rational  transfer  function  in  the  form 
//  Kg(s  +  zl)  (s  +  z2)  . . ./ (s+pl)  (s+p2)  .  .  . 

class  FT  { 

double  Kg;  //  high  frequency  gain. 

Polynom  numerator; 

Polynom  denominator; 

public : 

FT  0  {  } 

//  By  default: 

//  FT  (  const  FT  &  ); 

//  FT  &  operator  =  (  const  FT  &  )  ; 

int  getNumZeros  ()  const  {  return  numerator . getNumZeros () ;  } 

int  getNumZerosNotZero  ()  const  { 

return  numerator . getNumZerosNotZero {)  ;  } 

int  getNumPolos  ()  const  {  return  denominator . getNumZeros () ;  } 

int  getNumPolesNotZero  ()  const  { 

return  denominator . getNumZerosNotZero {) ;  } 

qcomplex  getZero  (int  i)  const  { 

return  numerator . getZero  ( i ) ;  } 

qcomplex  getZeroNotZero  (int  i)  const  { 

return  numerator . getZeroNotZero  ( i ) ;  } 

qcomplex  getPole  (int  i)  const  { 

return  denominator . getZero  ( i ) ;  } 

qcomplex  getPoleNotZero  (int  i)  const  { 

return  denominator . getZeroNotZero ( i ) ;  } 

void  addZero  (  qcomplex  z  ){  numerator . addZero ( z ) ;  } 
void  addPole  (  qcomplex  p  ){  denominator . addZero (p) ;  } 
void  setKg  (  double  k  ) {  Kg  =  k;  } 
double  getKg  ()  const  {  return  Kg;  } 
qcomplex  evaluatein  (  double  free  )  const  { 

return  ( std : : complex<double> )  numerator . evaluatein ( free)  *  Kg  / 
(std: : complex<double>)  denominator . evaluatein ( free) ;  } 

double  getArea  (double  wini,  double  wEnd)  const; 
friend  FT  operator  *  (  const  FT  &  ftl,  const  FT  &  ft2  ); 


FT  operator  *  (  const  FT  &  ftl,  const  FT  &  ft2  ); 


//  Represents  a  bound  at  a  certain  frequency: 

// 

//  It  is  a  double  set  of  values,  with  the  maximum  and  the  minimum 
//  permitted  values  for  each  phase  value. 

class  Bound  { 
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private : 


double 

double 

frequency; 
numValues ; 

// 

double 

*  minValues; 

// 

double 

*  maxValues; 

// 

phase  points  in  which  the  bound  has  been 
//  evaluated 

bounds  in  solid  line  in  QFT  of  MATLAB 
bounds  in  dashed  line  in  QFT  of  MATLAB 


//  NOTE:  there  is  no  relation  between  minValues  and  maxValues;  i.e., 
//  it  may  happen  minValues [ i ] >maxValues [ i ] 

double  getDeltaPhase ( )  const  {  return  360/ (numValues-1) ;  } 

//  gives  the  bound  at  a  certain  phase  (extrapolation  if  necessary) 
void  getBoundsForPhase (  double  phase,  double  &  min,  double  &  max  ) 
const ; 


public : 


Bound  0 {  minValues  =  maxValues  =  0; 
std: :vector 

Bound  (  std::istream  &  is  )  ; 

-Bound  ( ) ; 


//  necessary  to  use  later 

// 


Bound  (  const  Bound  &  ) ; 

Bound  &  operator  =  (  const  Bound  &  ) ; 

double  getFrequency ( )  const  {  return  frequency;  } 
double  getCost  (  const  FT  &  ftla  )  const; 

//  necessary  to  use  later  std::vector 

bool  operator<  (  const  Bound  &  b  )  const  {  return  true;  } 
bool  operator==  (  const  Bound  &  b  )  const  {  return  false;  } 


//  Set  of  bounds 
class  Bounds  { 

std::vector  <Bound>  Boundslist; 

Bounds  (  const  Bounds  &  ) ; 

Bounds  &  operator  =  (  const  Bounds  &  ) ; 
public : 

Bounds  (  std::istream  &  is  ); 

int  getNumBounds  ()  const  {  return  Boundslist . size () ;  } 

double  getCostDueToBound  (  int  i,  const  FT  &  ftla  )  const; 


//  Class  to  calculate  the  cost  of  a  controller.  This  cost  has  two 
//  components: 

//  -  basic  cost  (function  to  minimize) 
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//  -  penalty  (contraints  that  have  to  be  fulfilled) 

class  FCost  { 


double  aO,  al,  a2; 
double  wl,  w2 ; 

double  b; 

Bounds  *  theBounds; 
FT  thePlant; 


//  basic  cost  function  parameters 
//  interval  of  frequencies  to  measure  the 
//  controller  area 
//  penalty  parameter 
//  QFT  bounds 
/ /  QFT  nominal  plant 


FCost  (  const  FCost  &  ) ; 

FCost  &  operator  =  (  const  FCost  &  ) ; 


double  getBasicCost  (  const  FT  &  controller  )  const; 
double  getPenalty  (  const  FT  &  controller  )  const; 


public : 


FCost  0; 
-FCost  0; 


bool  fulfilBounds  (  const  FT  &  contr  )  const  { 
return  getPenalty (contr )  ==  0;  } 

double  getCost  (  const  FT  &  contr  )  const  { 

return  getBasicCost (contr)  +  getPenalty ( contr ) ;  } 


//  Utility  class 
class  XContr  { 

//  number  of  controller  fixed  poles,  i.e.,  the  algorithm 
/ /  does  not  move  them 
int  numFixedPoles ; 
double  *  theFixedPoles ; 

//  number  of  parameters  that  corresponds  to  the  controller 
//  zeros 
int  numZeros; 

//  the  parameters  must  be: 

//  -  the  first  is  the  controller  Kg 

//  -  the  following  "numZeros"  are  the  controller  zeros 

//  -  the  rest  are  the  controller  poles 

XContr  (  const  XContr  &  ) ; 

XContr  &  operator  =  (  const  XContr  &  ) ; 

public : 

XContr  ( ) ; 

-XContr  ( ) ; 

FT  getController  (  const  Phenotype  &  phenotype  ) ; 


#endif 
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//  BOUNDS. cpp  QFT  bounds  representation  and  controller  cost  calculation. 

// 

//  Authors:  Mario  Garcia-Sanz,  Carlos  Molins 
//  1 . September . 20 0 9 

// - 


#include  <math.h>  //  for  floor  () 

#include  <strstream> 

#include  <fstream> 

#include  "bounds. h" 


#define  M_PI  3.14159265358979323846 


//  Function  to  obtain  dB  and  phase 

void  dbphase  (  qcomplex  c,  double  &  db,  double  &  phase  ) { 
db=2  0*logl0  (  c.absO  ); 
phase=l 8  0  *c . arg ( ) /M_PI ; 

//  to  transform  from  (-180,180]  to  (-360,0] 
if  (phase  >  0)  phase  -=  360.0; 


void  Polynom: : addZero  (  qcomplex  c  ) { 

Zeroslist . push_back  (c) ; 

//  add  the  conjugate 

if  (  c  .  imag  ( )  )  Zeroslist .  push_back  (c.conjO); 

std::sort  (  Zeroslist . begin () ,  Zeroslist . end ( )  )  ; 

if  (  (c.realO  ==  0)  &&  (c.imagO  ==  0)  )  numOriginZeros+f; 


qcomplex  Polynom: : evaluatein  (  double  frequency  )  const  { 
std : : complex<double>  jw  (0,  frequency); 

std : : complex<double>  res  (1);  //  initialize  to  the  unity 

for  (int  i=0;  i<Zeroslist . size ( ) ;  i++) { 

res  *=  ( jw+ (std: : complex<double>) Zeroslist [i]  )  ; 

} 

return  res; 


Polynom  operator  *  (  const  Polynom  &  pi,  const  Polynom  &  p2  )  { 

Polynom  res; 

for  (int  i=0;  i<pl . Zeroslist . size () ;  i++)  res . Zeroslist . push_back 
(pl.Zeroslist[i] ) ; 

for  (int  i=0;  i<p2 . Zeroslist . size ()  ;  i++)  res . Zeroslist . push_back 
(p2 . Zeroslist [ i ] ) ; 
return  res; 


double  FT::getArea  (double  wini,  double  wEnd)  const  { 
//  area  of  ln|G(jw) | 
double  area  =  0; 
const  int  NUMW  =  100; 

double  deltaw  =  (wEnd  -  wIni)  /  NUMW; 
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double  w  =  wini; 

for  (int  1=0;  1<NUMW;  1++) { 

qcomplex  g  =  evaluatein  (w)  ; 
area  +=  deltaw  *  log  (g.absO); 
w  +=  deltaw; 

} 

return  area  >  0  ?  area  :  0.0; 

} 


FT  operator  *  (  const  FT  &  ftl,  const  FT  &  ft2  ) { 

FT  res; 

res.setKg  (  ftl. Kg  *  ft2.Kg  ); 

res . numerator  =  ftl . numerator  *  ft2 . numerator ; 

res . denominator  =  ftl . denominator  *  ft2 . denominator; 

return  res; 

} 


Bound: :Bound  (  std: :istream  &  is  ) { 
is  >>  frequency  >>  numValues; 
minValues  =  new  double  [numValues]; 
maxValues  =  new  double  [numValues]; 

for  (int  1=0;  i<numValues;  it+)  is  >>  minValues [ 1 ] ; 
for  (int  1=0;  i<numValues;  it+)  is  >>  maxValues [ 1 ] ; 

} 

Bound: : ~Bound  () { 

if  (minValues)  delete  []  minValues; 
if  (maxValues)  delete  []  maxValues; 

} 


Bound: :Bound  (  const  Bound  &  bnd  ) { 
frequency  =  bnd . frequency ; 
numValues  =  bnd . numValues ; 
minValues  =  new  double  [numValues]; 
maxValues  =  new  double  [numValues]; 
for  (int  1=0;  i<numValues;  i++)  { 
minValues [1]  =  bnd. minValues [1]  ; 
maxValues [1]  =  bnd. maxValues [1] ; 

} 


Bound  &  Bound :: operator  =  (  const  Bound  &  bnd  ) { 
if  (  &bnd  ! =  this  ) [ 

if  (minValues)  delete  []  minValues; 
if  (maxValues)  delete  []  maxValues; 
frequency  =  bnd . frequency ; 
numValues  =  bnd . numValues ; 
minValues  =  new  double  [numValues]; 
maxValues  =  new  double  [numValues]; 
for  (int  1=0;  i<numValues;  i++)  { 
minValues [i]  =  bnd . minValues [ 1 ] ; 
maxValues [i]  =  bnd . maxValues [ i ]  ; 

} 

} 

return  *this; 

} 

double  Bound :: getCost  (  const  FT  &  ftla  )  const  [ 
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qcomplex  c  =  f tla . evaluatein  {  getFrequency { )  ) ; 

double  db,  phase; 
dbphase  (  c,  db,  phase  ) ; 
double  min,  max; 

getBoundsForPhase (  phase,  min,  max  ); 
if  (  max  >  min  ) { 

if  (  db  <  min  )  return  min  -  db; 
if  (  db  >  max  )  return  db  -  max; 

} 

else  { 

if  (  (db  <  min)  &&  (db  >  max)  ) 

return  (min-db) > (db-max)  ?  db-max  :  min-db; 

} 

return  0.0; 


void  Bound :: getBoundsForPhase (  double  phase,  double  &  min, 

double  &  max  )  const  { 

//  linear  interpolation 

if  (  (phase>0)  | |  (phase<-360)  )  throw; 

double  delta  =  -  getDeltaPhase ( ) ; 

//  the  bounds  are  for  0 : -delta : -3 60  (MATLAB  notation) 
double  mudo  =  phase  /  delta; 
int  ind  =  (int)  floor  (mudo) ; 
if  (ind  ==  numValues-1)  ind — ; 

//  x-xl/x2-xl  =  Y-yl/y2-yl 

//  pointl  corresponds  to  ind,  point2  to  ind+1 

min  =  (  (minValues [ ind+1 ]  -  minValues [ind] )* (phase-ind*delta) /delta  ) 
+  minValues [ ind]  ; 

max  =  (  (maxValues [ ind+1 ]  -  maxValues [ind] )* (phase-ind*delta) /delta  ) 
+  maxValues [ind]; 


Bounds :: Bounds  (  std::istream  &  is  ) [ 
char  cadena  [30000]; 

while  (  !  is . eof ( ) )  { 

is . get line (cadena,  sizeof (cadena) ) ; 

std : : istrstream  iss  (cadena,  sizeof  (cadena)  ); 

Bound  bnd  (iss) ; 
if  (  iss . good ( )  )  { 

Boundslist . push_back  (bnd); 

} 

} 

} 

double  Bounds :: getCostDueToBound  (  int  i,  const  FT  &  ftla  )  const  { 
return  Boundslist [ i ]. getCost  (ftla); 


FCost: :FCost  () { 

std :: if stream  ifsConf  ("cfgcost.txt"); 
if  (  !  if sConf . good ( )  ){ 

std::cerr  <<  "\tCan  not  open  the  file  to  configure  the  fitness  function" 
<<  std : : endl ; 

throw; 

} 

char  cad  [250 ] ; 
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char  ficBounds  [250]; 
ifsConf  >>  ficBounds; 
if sConf . getline  (cad,  sizeof (cad) ) ; 


ifsConf  >>  aO; 
if sConf . getline 
ifsConf  >>  al; 
ifsConf . getline 
ifsConf  >>  a2; 
ifsConf . getline 
ifsConf  >>  wl; 
ifsConf . getline 
ifsConf  >>  w2; 
ifsConf . getline 
ifsConf  >>  b; 
ifsConf . getline 


(cad,  sizeof (cad) ) ; 
(cad,  sizeof (cad) ) ; 
(cad,  sizeof (cad) ) ; 
(cad,  sizeof (cad) ) ; 
(cad,  sizeof (cad) ) ; 
(cad,  sizeof (cad) ) ; 


double  k; 
ifsConf  >>  k; 

if sConf . getline  (cad,  sizeof (cad) ) ; 
thePlant . setKg  (  k  ); 


int  numZeros; 
ifsConf  >>  numZeros; 
if sConf . getline  (cad,  sizeof (cad) ) ; 
for  (int  i  =  0;  i  <  numZeros;  i+t) { 
qcomplex  zero; 
ifsConf  >>  zero; 

if sConf . getline  (cad,  sizeof (cad) ) ; 
thePlant . addZero  (  zero  ); 

} 


int  numPoles; 
ifsConf  >>  numPoles; 
if sConf . getline  (cad,  sizeof (cad) ) ; 
for  (int  i  =  0;  i  <  numPoles;  itt) { 
qcomplex  pole; 
ifsConf  >>  pole; 

if sConf . getline  (cad,  sizeof (cad) ) ; 
thePlant . addPole  (  pole  ); 

} 


std :: if stream  ifs  (ficBounds); 
if  (  !  ifs . good ( )  ) { 

std::cerr  <<  "\tCan  not  open  the  bounds  file"  <<  std::endl; 
throw; 


theBounds  =  new  Bounds  (ifs) ; 

} 


FCost: :~FCost  () { 
delete  theBounds; 

} 

double  FCost :: getBasicCost  (  const  FT  &  controller  )  const  { 

//  Fitness  function  based  on  Thompson  and  Nwokah[38]  and  Thompson  [39] 
double  costO  =  aO  *  controller . getKg ( )  *  controller . getKg () ; 

//  for  the  costl  remove  the  less  significative  poles 
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//  until  np=nz 
double  c  =  0; 

//  costl  takes  into  account  poles  and  zeros  different  from  0 

int  numZeros  =  controller . getNumZerosNotZero {) ; 
if  (controller . getNumPolesNotZero { )  <  numZeros)  numZeros  = 
controller . getNumPolesNotZero { ) ; 
for  (int  1=0;  1  <  numZeros;  i++) { 

//  It  is  supossed  that  the  controller  has  no  complex  zeros  and  no 
//  complex  poles. 

qcomplex  pi  =  controller . getPoleNotZero ( 1 ) ; 
qcomplex  zi  =  controller . getZeroNotZero ( 1 ) ; 
double  dif  =  pi . real () -zi . real () ; 

c  +=  (  (dif  *  dif  *  dif  *  dif)  +  1  )  /  (pi. real  ()  *  zi.realO); 

} 

double  costl  =  al  *  c; 

double  cost2  =  a2  *  controller . getArea (wl ,  w2); 
return  costO  +  costl  +  cost2; 


double  FCost : : getPenalty  (  const  FT  &  controller  )  const  { 
double  cost  =  0; 

FT  ftla  =  controller  *  thePlant; 

int  numBounds  =  theBounds->getNumBounds  (); 

for  (int  1=0;  i<numBounds;  i++)  { 

double  c  =  theBounds->getCostDueToBound  (  1,  ftla  ); 
cost  +=  b  *  c  *  c; 

} 

return  cost; 

} 


XContr: :XContr  () { 

std :: if stream  ifsConf  ("cfgcontr.txt"); 
if  (  !  if sConf . good ( )  ){ 

std::cerr  <<  "\tCan  not  open  the  controller  configuration  file" 
<<  std : : endl ; 

throw; 

} 

char  cad  [250 ] ; 

ifsConf  >>  numFixedPoles; 
if sConf . getline  (cad,  sizeof (cad) ) ; 
theFixedPoles  =  new  double  [numFixedPoles]; 
for  (int  1=0;  1  <  numFixedPoles;  i++) { 
ifsConf  >>  theFixedPoles[i]; 
if sConf . getline  (cad,  sizeof (cad) ) ; 

} 

ifsConf  >>  numZeros; 

if sConf . getline  (cad,  sizeof (cad)); 


XContr: :~XContr  () { 

delete  []  theFixedPoles; 

} 

FT  XContr :: getController  (  const  Phenotype  &  phenotype  ){ 
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//  NOTE:  the  number  of  poles  different  from  0  to  optimize  hve  to  be 
//  the  same  as  the  number  of  zeros  different  from  0  to  optimize. The 
//  remainder  poles  or  zeros  have  to  be  fixed  (Thompson) 

FT  controller; 

controller . setKg  (  phenotype . getParameterValue  (0)  ); 

for  (int  i  =  1;  i  <=  numZeros;  itt) { 

controller . addZero  (  (qcomplex)  phenotype . getParameterValue  (1)  ); 

} 

int  numParameters  =  Parameters :: getNumParameters {) ; 
for  (int  i  =  numZerostl;  i  <  numParameters;  i++) { 

controller . addPole  (  (qcomplex)  phenotype . getParameterValue  (1)  ); 

} 

for  (int  1=0;  1  <  numFixedPoles ;  i++) { 

controller . addPole  (  (qcomplex)  theFixedPoles [ 1 ]  ); 

} 

return  controller; 


static  FCost  fcost; 
static  XContr  xcontr; 


//  Fitness  function 

double  externalEvaluation  (  const  Phenotype  &  phenotype  ) { 

FT  controller  =  xcontr . getController  (phenotype); 

//  the  objective  of  the  Genetic  Algorithms  is  to  find  the  maximum 
return  1/fcost . getCost (  controller  ); 


//  Function  to  check  the  feasibility 
//  (no  constraint  violation) 

bool  isValidExternalSolution  (  const  Phenotype  &  phenotype  ) { 
//  if  no  constraint: 

//  return  true; 

//  If  there  are  constraints: 

FT  controller  =  xcontr . getController  (phenotype); 
return  f cost . f ulf ilBounds (  controller  ); 


//  OPTQFT.cpp  This  is  the  main  program  that  tries  to  obtain  an  optimal 

//  controller  applying  the  Genetic  Algorithms. 

// 

//  Authors:  Mario  Garcia-Sanz,  Carlos  Molins 
//  1 . September . 20 0 9 

// - 

#lnclude  <math.h> 

#lnclude  <iostream> 

#lnclude  <fstream> 

#lnclude  "aagg.h" 


using  namespace  std; 
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void  main (  int  argc,  char  *  argv  []  ) { 

//  Read  configuration 
ifstream  ifsConf  ("cfg.txt"); 
if  (  !  if sConf . good ( )  ){ 

cerr  <<  "\tCan  not  open  the  configuration  file"  <<  endl; 
throw; 

} 

char  cad  [250 ] ; 

int  numParameters ; 

ifsConf  >>  numParameters; 

if sConf . getline  (cad,  sizeof (cad) ) ; 

double  *  minParam  =  new  double  [numParameters]; 
for  (int  i  =  0;  i  <  numParameters;  i++) { 
ifsConf  >>  minParam[i]; 
if sConf . getline  (cad,  sizeof (cad) ) ; 

} 

double  *  maxParam  =  new  double  [numParameters]; 
for  (int  i  =  0;  i  <  numParameters;  i++) [ 
ifsConf  >>  maxParam[i]; 
if sConf . getline  (cad,  sizeof (cad) ) ; 

} 

double  precision; 

ifsConf  >>  precision; 

if sConf . getline  (cad,  sizeof (cad) ) ; 

int  numindividuals ; 
double  pCross,  pMutate; 
int  showContinuous ; 
int  maxGenerations ; 
double  minFitness; 

ifsConf  >>  numindividuals; 

if sConf . getline  (cad,  sizeof (cad) ) ; 

ifsConf  >>  pCross; 

if sConf . getline  (cad,  sizeof (cad) ) ; 
ifsConf  >>  pMutate; 

if sConf . getline  (cad,  sizeof (cad) ) ; 

ifsConf  >>  showContinuous; 

if sConf . getline  (cad,  sizeof (cad) ) ; 

ifsConf  >>  maxGenerations; 

if sConf . getline  (cad,  sizeof (cad) ) ; 

ifsConf  >>  minFitness; 

if sConf . getline  (cad,  sizeof (cad) ) ; 

ofstream  ofs  ("aagg.log"); 
if  (  !  of s . good ( )  )  { 

cerr  <<  "\tCan  not  open  the  output  file"  <<  endl; 
return ; 

} 

//  Initialize  the  parameters:  numGenes,  valueMin,  valueMax 

//  precision: 

// 
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//  The  precision  is  calculated  by: 

//  r=2 * *numGenes  -1 
//  prec= (max-min) /r 

for  (int  i  =  0;  i  <  numParameters ;  i++) { 

double  y  =  ( (maxParam [ i ] -minParam [ i ]) /precision) +1  ; 

int  numGenes  =  floor  (  log ( (double) y) /log ( (double) 2 )  )  +  1; 

Parameters :: addParameter  (  numGenes,  minParam [1],  maxParam [1]  ); 

} 

//  Remove  the  file  to  obtain  evolution  graphics 
std :: of stream  ofsPriv  ("privl.txt"); 
ofsPr iv. close  0 ; 

//  Initialize  the  algorithm 
GeneticAlgorithm  ag; 

ag . initialize  (  numindividuals ,  pCross,  pMutate,  (bool)  showContinuous , 
maxGenerations ,  minFitness  )  ; 

//  execute 

ag .mainAlgorithm  (  ofs  )  ; 

delete  []  minParam; 
delete  []  maxParam; 


//  RANDOM. h  Functions  with  random  numbers. C++  adaptation  of  the 

//  Goldberg  basic  code  [23]  . 

// 

// 

//  Authors:  Mario  Garcia-Sanz,  Carlos  Molins 
//  1 . September . 20 0 9 

// - 


#if  [defined  _ RANDOM_H _ 

#define  _ RANDOM_H _ 

/ /  random  number  between  1  and  0 
double  randOl  (); 

//  gives  TRUE  with  probability  p 
bool  flip  (  double  p) ; 

//  gives  a  random  index  (integer)  between  min  and  max 
int  getRandomIndex  (int  min,  int  max) ; 

#endif 


// 

// 

// 

// 

// 

// 

// 


RANDOM. cpp  Functions  with  random  numbers. C++ 

Goldberg  basic  code  [23] . 


Authors:  Mario  Garcia-Sanz,  Carlos  Molins 
1 . September .2009 


adaptation 


of  the 
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#include  <stdlib.h> 
#include  "random. h" 


double  randOl  (){ 

/ /  random  number  between  0  and  1 
int  r  =  rand ( ) ; 

return  (double)  r  /  (double)  RAND_MAX; 

} 


bool  flip  (double  p) { 

//  gives  TRUE  with  probability  p 
return  randOlO  <=  p; 

} 


int  getRandomIndex  (int  min,  int  max) { 

//  gives  a  random  index  (integer)  between  min  and  max 

int  range  =  max  -  min  +  1  ; 

int  newMax  =  (long)  RAND_MAX  +  1  -  (  (long)  RAND_MAX  +  1)  %  range; 

int  temp  =  rand(); 

while  (temp  >  newMax  )  temp  =  rand ( ) ; 
return  temp%range  +  min; 
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8.3.  TEXT  FILES 


CFG.txt  This  file  is  to  define  range  of  variation  of  controller 

parameters,  as  well  as  the  crossover  and  mutation  probabilities, 
maximum  number  of  generations,  the  number  of  individuals  in  each 
generation, ... 


11  number  of  parameters 

0.01  minimum  value  for  the  parameter  1 

0.01  minimum  value  for  the  parameter  2 

0.01  minimum  value  for  the  parameter  3 

0.01  minimum  value  for  the  parameter  4 

0.01  minimum  value  for  the  parameter  5 

0.01  minimum  value  for  the  parameter  6 


0 .01 

minimum 

value 

for 

the 

parameter 

7 

0 .01 

minimum 

value 

for 

the 

parameter 

8 

0 .01 

minimum 

value 

for 

the 

parameter 

9 

0 .01 

minimum 

value 

for 

the 

parameter 

10 

0 .01 

minimum 

value 

for 

the 

parameter 

11 

1000 

maximum 

value 

for 

the 

parameter 

1 

5 

maximum 

value 

for 

the 

parameter 

2 

5 

maximum 

value 

for 

the 

parameter 

3 

5 

maximum 

value 

for 

the 

parameter 

4 

10 

maximum 

value 

for 

the 

parameter 

5 

10 

maximum 

value 

for 

the 

parameter 

6 

1000 

maximum 

value 

for 

the 

parameter 

7 

1000 

maximum 

value 

for 

the 

parameter 

8 

1000 

maximum 

value 

for 

the 

parameter 

9 

1000 

maximum 

value 

for 

the 

parameter 

10 

1000 

maximum 

value 

for 

the 

parameter 

11 

.01 

parameters  precision 

100  number  of  individuals 

0.1  crossover  probability 

0.01  mutation  probability 

1  show  continuously  (1)  or  only  (0) 

500  maximum  number  of  generations 

100000  minimum  desired  fitness 


CFGCONTR.txt 


This  file  is  to  define  some  controller  characteristics. 


0 

5 


number  of  controller  fixed  poles 

number  of  parameters  that  corresponds  to  controller  zeros 


CFGCOST . txt 


This  file  is  to  define  the  fitness  function  parameters 
values,  as  well  as  the  nominal  plant. 


bnd_n 
0 .001 
0 .00001 
1 

0 .01 
1 
1 

1440 


bounds  file 
aO  parameter 
al  parameter 
a2  parameter 
wlni  to  measure 
wEnd  to  measure 
a3  parameter 
gain  of  the  QFT 


the  "excess  controller  bandwidth" 
the  "excess  controller  bandwidth" 

nominal  plant 
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1 

0.5 

3 

0 

10 

(4. 8, 3. 6) 


number  of  zeros  of  the  QFT  nominal  plant 

nominal  plant  zero;  for  complex:  (2,3)  only  one  of  the  conjugates 

number  of  poles  of  the  QFT  nominal  plant  (without  conjugates) 
nominal  plant  pole;  for  complex:  (2,3)  only  one  of  the  conjugates 

nominal  plant  pole;  for  complex:  (2,3)  only  one  of  the  conjugates 

nominal  plant  pole;  for  complex:  (2,3)  only  one  of  the  conjugates 


8.4.  HOW  TO  RUN 

First  of  all,  move  all  the  .cpp  and  .h  files  to  the  folder: 

C:\Program  FilesVMicrosoft  Visual  Studio  9.0\VC 

Once  this  is  done,  write  at  the  Visual  Studio  2008  Command  Prompt: 

cl  -EHsc  OPTQFT.cpp  aagg.cpp  random.cpp  bounds. cpp 

With  this  you  get:  optqft.exe 

Then,  move  this  optqft.exe  to  the  folder  where  the  Matlab  and  text  files  are. 

Finally,  read  the  files  just  in  case  any  other  configuration  is  desired,  and  run  at  the  Matlab 
Command  Window:  Autoloopshaping_AG.m 
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